options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(1,4))
      drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
      drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
      drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
      drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}

regression

ex5.stan

multi normal regression

data{
  int N;
  int K;
  vector[N] y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
  real<lower=0> s;
}
model{
  vector[N] m=X*b;
  y~normal(m,s);
}
generated quantities{
  vector[N] y1;
  vector[N] m1=X*b;
  for(i in 1:N){
    y1[i]=normal_rng(m1[i],s);
  }
}

normal linear models

n=30
mdl=cmdstan_model('./ex5.stan')

single regression

tb=tibble(x=runif(n,0,9),
          y=rnorm(n,x,1))
f=formula(y~x)
par(mfrow=c(1,1))
plot(tb$x,tb$y)

qplot(data=tb,x,y)

X=model.matrix(f,tb)

data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
## Initial log joint probability = -798.475 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       16      -14.6353   6.65847e-05   0.000885314      0.8846      0.8846       22    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
mle
##  variable estimate
##    lp__     -14.64
##    b[1]       0.83
##    b[2]       0.82
##    s          0.99
##    y1[1]      5.89
##    y1[2]      4.33
##    y1[3]      2.59
##    y1[4]      2.39
##    y1[5]      3.49
##    y1[6]      5.18
##    y1[7]      6.10
##    y1[8]      2.58
##    y1[9]      7.74
##    y1[10]     3.37
##    y1[11]     2.70
##    y1[12]     2.57
##    y1[13]     3.50
##    y1[14]     0.49
##    y1[15]     7.64
##    y1[16]     7.02
##    y1[17]     0.22
##    y1[18]     6.68
##    y1[19]     4.37
##    y1[20]     5.99
##    y1[21]     1.46
##    y1[22]     7.57
##    y1[23]     5.39
##    y1[24]     6.89
##    y1[25]     6.96
##    y1[26]     4.34
##    y1[27]     3.90
##    y1[28]     8.82
##    y1[29]     3.70
##    y1[30]     6.12
##    m1[1]      5.57
##    m1[2]      4.32
##    m1[3]      1.75
##    m1[4]      2.84
##    m1[5]      4.06
##    m1[6]      4.33
##    m1[7]      5.36
##    m1[8]      1.68
##    m1[9]      6.99
##    m1[10]     3.06
##    m1[11]     2.84
##    m1[12]     3.55
##    m1[13]     4.85
##    m1[14]     1.74
##    m1[15]     7.21
##    m1[16]     7.95
##    m1[17]     1.72
##    m1[18]     8.00
##    m1[19]     2.53
##    m1[20]     6.36
##    m1[21]     3.30
##    m1[22]     7.04
##    m1[23]     6.98
##    m1[24]     6.44
##    m1[25]     6.08
##    m1[26]     4.66
##    m1[27]     5.35
##    m1[28]     8.02
##    m1[29]     3.37
##    m1[30]     6.67
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -16.17 -15.84 1.23 1.04 -18.63 -14.84 1.01      539     1012
##    b[1]     0.83   0.83 0.42 0.41   0.13   1.51 1.01      520      527
##    b[2]     0.82   0.82 0.08 0.08   0.69   0.95 1.00      623      773
##    s        1.07   1.05 0.15 0.14   0.86   1.34 1.00     1290     1176
##    y1[1]    5.56   5.55 1.07 1.05   3.84   7.30 1.00     2074     1969
##    y1[2]    4.31   4.30 1.08 1.04   2.54   6.06 1.00     1853     1931
##    y1[3]    1.75   1.75 1.13 1.09  -0.19   3.60 1.00     1814     1781
##    y1[4]    2.80   2.79 1.11 1.09   0.95   4.63 1.00     2038     1768
##    y1[5]    4.07   4.06 1.11 1.08   2.25   5.86 1.00     1808     1894
##    y1[6]    4.32   4.33 1.11 1.03   2.45   6.13 1.00     1831     1828
##    y1[7]    5.35   5.35 1.10 1.08   3.54   7.16 1.00     1924     1996
##    y1[8]    1.72   1.72 1.12 1.12  -0.16   3.52 1.00     1463     1757
##    y1[9]    7.00   6.98 1.09 1.07   5.23   8.77 1.00     2117     2015
##    y1[10]   3.08   3.08 1.13 1.05   1.22   4.93 1.00     1945     1866
##    y1[11]   2.87   2.86 1.11 1.07   1.07   4.74 1.00     1947     1846
##    y1[12]   3.56   3.54 1.11 1.07   1.73   5.40 1.00     1900     1945
##    y1[13]   4.87   4.87 1.10 1.11   3.07   6.63 1.00     2030     1731
##    y1[14]   1.74   1.71 1.13 1.13  -0.08   3.58 1.00     1558     1769
##    y1[15]   7.21   7.20 1.10 1.12   5.41   9.04 1.00     1958     1929
##    y1[16]   7.96   7.96 1.16 1.13   6.06   9.84 1.00     1607     1601
##    y1[17]   1.72   1.71 1.14 1.11  -0.10   3.59 1.00     1646     1840
##    y1[18]   7.99   7.99 1.13 1.09   6.13   9.89 1.00     2186     1881
##    y1[19]   2.50   2.49 1.10 1.07   0.66   4.34 1.00     1537     1897
##    y1[20]   6.36   6.38 1.10 1.07   4.50   8.13 1.00     2212     1859
##    y1[21]   3.27   3.30 1.12 1.06   1.43   5.08 1.00     1955     1721
##    y1[22]   7.03   7.03 1.14 1.17   5.22   8.93 1.00     1940     1772
##    y1[23]   6.99   6.98 1.13 1.10   5.18   8.78 1.00     2043     1858
##    y1[24]   6.42   6.40 1.13 1.09   4.59   8.30 1.00     1751     1971
##    y1[25]   6.14   6.13 1.12 1.10   4.33   7.97 1.00     2187     1930
##    y1[26]   4.66   4.65 1.10 1.10   2.85   6.50 1.00     1917     1770
##    y1[27]   5.36   5.39 1.11 1.10   3.53   7.10 1.00     1900     1772
##    y1[28]   8.02   8.00 1.15 1.11   6.14   9.89 1.00     1962     1951
##    y1[29]   3.40   3.41 1.09 1.05   1.67   5.19 1.00     1900     1828
##    y1[30]   6.70   6.70 1.08 1.04   4.92   8.42 1.00     1987     1815
##    m1[1]    5.57   5.57 0.21 0.21   5.23   5.92 1.00     2126     1452
##    m1[2]    4.32   4.32 0.20 0.20   3.99   4.65 1.00      953     1037
##    m1[3]    1.76   1.76 0.35 0.34   1.19   2.31 1.01      527      534
##    m1[4]    2.84   2.84 0.27 0.27   2.40   3.27 1.01      570      705
##    m1[5]    4.07   4.07 0.21 0.20   3.73   4.41 1.01      824      941
##    m1[6]    4.33   4.34 0.20 0.20   4.01   4.67 1.00      963     1123
##    m1[7]    5.37   5.37 0.20 0.20   5.03   5.70 1.00     2033     1474
##    m1[8]    1.68   1.69 0.35 0.35   1.11   2.25 1.01      526      534
##    m1[9]    6.99   6.99 0.29 0.28   6.52   7.47 1.00     1758     1408
##    m1[10]   3.06   3.06 0.25 0.26   2.64   3.47 1.01      590      770
##    m1[11]   2.84   2.84 0.27 0.27   2.40   3.28 1.01      570      705
##    m1[12]   3.55   3.56 0.23 0.23   3.17   3.93 1.01      665      915
##    m1[13]   4.85   4.85 0.20 0.20   4.53   5.18 1.00     1428     1397
##    m1[14]   1.74   1.74 0.35 0.34   1.17   2.30 1.01      527      531
##    m1[15]   7.21   7.21 0.30 0.30   6.73   7.72 1.00     1630     1398
##    m1[16]   7.96   7.95 0.36 0.36   7.37   8.55 1.00     1300     1365
##    m1[17]   1.73   1.73 0.35 0.34   1.16   2.29 1.01      527      534
##    m1[18]   8.00   8.00 0.36 0.36   7.41   8.60 1.00     1287     1365
##    m1[19]   2.53   2.53 0.29 0.29   2.05   2.99 1.01      551      641
##    m1[20]   6.36   6.36 0.25 0.24   5.96   6.77 1.01     2038     1428
##    m1[21]   3.30   3.30 0.24 0.24   2.90   3.69 1.01      621      893
##    m1[22]   7.04   7.04 0.29 0.29   6.57   7.53 1.00     1730     1388
##    m1[23]   6.98   6.98 0.29 0.28   6.51   7.46 1.00     1765     1408
##    m1[24]   6.45   6.45 0.25 0.24   6.04   6.87 1.01     2010     1427
##    m1[25]   6.08   6.08 0.23 0.23   5.70   6.47 1.01     2117     1400
##    m1[26]   4.66   4.66 0.20 0.20   4.34   4.99 1.00     1219     1313
##    m1[27]   5.35   5.35 0.20 0.20   5.01   5.68 1.00     2022     1473
##    m1[28]   8.02   8.02 0.36 0.36   7.43   8.62 1.00     1282     1365
##    m1[29]   3.37   3.37 0.24 0.24   2.97   3.76 1.01      632      901
##    m1[30]   6.67   6.67 0.27 0.26   6.24   7.12 1.00     1920     1348

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)

qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)

par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(1:n,y,data=tb,col=I('red'))+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

multiple regression

tb=tibble(x1=runif(n,0,9),x2=runif(n,0,9),
          y=rnorm(n,x1-x2,1))
f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
## Initial log joint probability = -111.072 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       31       -9.1604    1.6178e-05   0.000692474           1           1       42    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__      -9.16
##    b[1]      -0.18
##    b[2]       1.06
##    b[3]      -0.99
##    s          0.82
##    y1[1]      2.99
##    y1[2]     -1.19
##    y1[3]      4.02
##    y1[4]     -5.16
##    y1[5]      5.39
##    y1[6]     -1.92
##    y1[7]     -0.87
##    y1[8]     -2.35
##    y1[9]     -1.63
##    y1[10]    -3.67
##    y1[11]    -7.58
##    y1[12]    -0.71
##    y1[13]    -0.74
##    y1[14]    -0.21
##    y1[15]    -2.75
##    y1[16]    -6.00
##    y1[17]    -0.03
##    y1[18]     0.40
##    y1[19]     4.35
##    y1[20]     0.18
##    y1[21]    -6.89
##    y1[22]    -2.82
##    y1[23]     2.08
##    y1[24]     0.14
##    y1[25]    -5.08
##    y1[26]    -4.82
##    y1[27]    -2.77
##    y1[28]     0.38
##    y1[29]     5.09
##    y1[30]     4.05
##    m1[1]      3.22
##    m1[2]     -0.57
##    m1[3]      4.13
##    m1[4]     -4.61
##    m1[5]      5.12
##    m1[6]     -0.99
##    m1[7]     -0.08
##    m1[8]     -2.43
##    m1[9]     -1.14
##    m1[10]    -3.57
##    m1[11]    -6.95
##    m1[12]    -1.69
##    m1[13]    -0.78
##    m1[14]     0.33
##    m1[15]    -3.46
##    m1[16]    -7.04
##    m1[17]     0.05
##    m1[18]    -0.70
##    m1[19]     3.77
##    m1[20]    -0.60
##    m1[21]    -6.09
##    m1[22]    -3.89
##    m1[23]     1.19
##    m1[24]     0.43
##    m1[25]    -5.64
##    m1[26]    -4.16
##    m1[27]    -0.64
##    m1[28]    -1.17
##    m1[29]     5.12
##    m1[30]     3.67
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median   sd  mad     q5   q95 rhat ess_bulk ess_tail
##    lp__   -11.38 -11.08 1.39 1.22 -14.05 -9.72 1.00      738     1115
##    b[1]    -0.20  -0.20 0.42 0.41  -0.89  0.48 1.01      382      413
##    b[2]     1.06   1.06 0.06 0.06   0.95  1.16 1.00      815     1217
##    b[3]    -0.99  -0.99 0.06 0.06  -1.09 -0.88 1.00      682     1101
##    s        0.91   0.90 0.12 0.12   0.73  1.14 1.00     1021      964
##    y1[1]    3.25   3.24 0.96 0.93   1.68  4.82 1.00     2206     2039
##    y1[2]   -0.56  -0.57 0.95 0.95  -2.11  1.01 1.00     1706     1845
##    y1[3]    4.12   4.13 0.96 0.93   2.54  5.66 1.00     1922     1745
##    y1[4]   -4.59  -4.60 0.96 0.94  -6.13 -3.01 1.00     1677     1599
##    y1[5]    5.12   5.13 0.98 0.96   3.52  6.70 1.00     1760     1783
##    y1[6]   -1.03  -1.03 1.01 0.99  -2.67  0.65 1.00     1260     1688
##    y1[7]   -0.08  -0.07 0.98 0.99  -1.71  1.47 1.00     1705     1820
##    y1[8]   -2.40  -2.43 0.92 0.83  -3.86 -0.86 1.00     1880     1746
##    y1[9]   -1.13  -1.11 0.95 0.88  -2.71  0.43 1.00     1963     1842
##    y1[10]  -3.57  -3.57 0.96 0.92  -5.19 -2.04 1.00     1967     1994
##    y1[11]  -6.95  -6.97 0.97 0.93  -8.50 -5.33 1.00     2087     1769
##    y1[12]  -1.69  -1.66 1.00 0.95  -3.42 -0.10 1.00     1719     2016
##    y1[13]  -0.80  -0.81 0.96 0.95  -2.43  0.75 1.00     2110     1831
##    y1[14]   0.30   0.28 0.97 0.95  -1.23  1.94 1.00     1858     1973
##    y1[15]  -3.46  -3.49 0.94 0.91  -4.96 -1.88 1.00     1890     1893
##    y1[16]  -7.04  -7.03 0.97 0.96  -8.61 -5.42 1.00     1693     1810
##    y1[17]   0.07   0.05 0.93 0.93  -1.45  1.51 1.00     1935     1673
##    y1[18]  -0.74  -0.73 0.93 0.92  -2.20  0.78 1.00     1961     1774
##    y1[19]   3.76   3.80 0.97 0.95   2.11  5.31 1.00     1909     1819
##    y1[20]  -0.60  -0.58 0.95 0.92  -2.15  0.97 1.00     1812     1892
##    y1[21]  -6.08  -6.07 0.98 0.99  -7.70 -4.53 1.00     1946     1972
##    y1[22]  -3.90  -3.89 0.97 0.97  -5.54 -2.33 1.00     1982     1855
##    y1[23]   1.17   1.14 0.94 0.92  -0.40  2.69 1.00     1513     1874
##    y1[24]   0.38   0.38 1.01 1.00  -1.28  2.11 1.00     1299     1711
##    y1[25]  -5.65  -5.63 0.96 0.93  -7.20 -4.09 1.00     1911     1951
##    y1[26]  -4.14  -4.11 0.96 0.93  -5.74 -2.58 1.00     1975     1659
##    y1[27]  -0.62  -0.61 0.97 0.96  -2.26  0.89 1.00     1684     2033
##    y1[28]  -1.19  -1.19 0.92 0.90  -2.69  0.28 1.00     1915     1854
##    y1[29]   5.09   5.09 1.01 0.96   3.47  6.78 1.00     1952     2042
##    y1[30]   3.69   3.70 0.95 0.93   2.11  5.22 1.00     1724     1802
##    m1[1]    3.23   3.24 0.29 0.28   2.75  3.71 1.00     1819     1619
##    m1[2]   -0.58  -0.58 0.28 0.28  -1.05 -0.11 1.01      401      520
##    m1[3]    4.14   4.15 0.33 0.32   3.59  4.68 1.00     1760     1563
##    m1[4]   -4.62  -4.62 0.27 0.26  -5.07 -4.18 1.01      700     1159
##    m1[5]    5.12   5.11 0.35 0.34   4.56  5.69 1.00     1199     1401
##    m1[6]   -1.01  -1.01 0.35 0.34  -1.58 -0.43 1.01      385      445
##    m1[7]   -0.09  -0.09 0.27 0.27  -0.52  0.35 1.01      418      547
##    m1[8]   -2.44  -2.44 0.19 0.20  -2.77 -2.14 1.01      722     1069
##    m1[9]   -1.13  -1.13 0.22 0.22  -1.50 -0.76 1.00     1705     1385
##    m1[10]  -3.57  -3.57 0.24 0.23  -3.95 -3.16 1.00     1941     1740
##    m1[11]  -6.95  -6.94 0.34 0.33  -7.50 -6.40 1.00     1822     1695
##    m1[12]  -1.68  -1.68 0.31 0.30  -2.18 -1.17 1.00      870     1286
##    m1[13]  -0.78  -0.77 0.24 0.24  -1.18 -0.37 1.00     1280     1347
##    m1[14]   0.33   0.33 0.22 0.21  -0.03  0.70 1.00     1840     1596
##    m1[15]  -3.45  -3.46 0.25 0.24  -3.85 -3.03 1.00     1651     1789
##    m1[16]  -7.04  -7.03 0.33 0.33  -7.59 -6.50 1.00     1742     1511
##    m1[17]   0.06   0.07 0.26 0.26  -0.38  0.51 1.00     1151     1595
##    m1[18]  -0.70  -0.69 0.17 0.16  -0.97 -0.42 1.00     1693     1420
##    m1[19]   3.78   3.79 0.34 0.33   3.21  4.33 1.00     1426     1545
##    m1[20]  -0.61  -0.61 0.24 0.25  -1.01 -0.21 1.01      432      662
##    m1[21]  -6.10  -6.09 0.31 0.29  -6.61 -5.61 1.00     1042     1255
##    m1[22]  -3.89  -3.89 0.22 0.22  -4.24 -3.53 1.00     1968     1639
##    m1[23]   1.18   1.19 0.28 0.28   0.73  1.63 1.01      461      718
##    m1[24]   0.42   0.42 0.38 0.37  -0.21  1.04 1.01      389      394
##    m1[25]  -5.64  -5.64 0.28 0.27  -6.10 -5.19 1.00     1748     1438
##    m1[26]  -4.16  -4.17 0.22 0.22  -4.53 -3.80 1.00     1437     1484
##    m1[27]  -0.63  -0.62 0.34 0.33  -1.21 -0.07 1.00      772     1253
##    m1[28]  -1.17  -1.17 0.17 0.17  -1.46 -0.90 1.00      953     1197
##    m1[29]   5.12   5.12 0.35 0.35   4.56  5.70 1.00     1664     1420
##    m1[30]   3.67   3.66 0.29 0.29   3.18  4.15 1.00      926     1380

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)

qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(1:n,y,data=tb,col=I('red'))+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

ANOVA

tb=tibble(c=sample(c('a','b','c'),n,replace=T),
          y=rnorm(n,(c=='b')*2-(c=='c')*2,1))
f=formula(y~c)
par(mfrow=c(1,1))
boxplot(y~c,tb)

qplot(data=tb,c,y,geom='boxplot')

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
## Initial log joint probability = -311.803 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       19      -12.4643   5.40671e-05   0.000515921      0.9009      0.9009       22    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -12.46
##    b[1]       0.15
##    b[2]       2.05
##    b[3]      -2.02
##    s          0.92
##    y1[1]      3.55
##    y1[2]     -1.81
##    y1[3]      0.17
##    y1[4]      0.73
##    y1[5]      1.75
##    y1[6]      1.00
##    y1[7]     -1.14
##    y1[8]      0.64
##    y1[9]     -0.77
##    y1[10]     1.05
##    y1[11]     0.94
##    y1[12]     1.39
##    y1[13]     0.22
##    y1[14]     0.62
##    y1[15]     1.26
##    y1[16]     1.60
##    y1[17]     1.74
##    y1[18]     2.78
##    y1[19]    -0.82
##    y1[20]     1.84
##    y1[21]    -2.87
##    y1[22]     1.70
##    y1[23]    -2.17
##    y1[24]    -2.73
##    y1[25]     1.83
##    y1[26]    -1.98
##    y1[27]     2.11
##    y1[28]     1.48
##    y1[29]     2.02
##    y1[30]    -0.82
##    m1[1]      2.19
##    m1[2]     -1.87
##    m1[3]      0.15
##    m1[4]      2.19
##    m1[5]      0.15
##    m1[6]      2.19
##    m1[7]      0.15
##    m1[8]      0.15
##    m1[9]      0.15
##    m1[10]     0.15
##    m1[11]     0.15
##    m1[12]     0.15
##    m1[13]     0.15
##    m1[14]     0.15
##    m1[15]     0.15
##    m1[16]     2.19
##    m1[17]     2.19
##    m1[18]     2.19
##    m1[19]     0.15
##    m1[20]     2.19
##    m1[21]    -1.87
##    m1[22]     2.19
##    m1[23]    -1.87
##    m1[24]    -1.87
##    m1[25]     2.19
##    m1[26]    -1.87
##    m1[27]     2.19
##    m1[28]     2.19
##    m1[29]     2.19
##    m1[30]    -1.87
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -14.78 -14.46 1.52 1.33 -17.79 -12.95 1.00      736      992
##    b[1]     0.16   0.17 0.29 0.27  -0.32   0.65 1.00      980     1057
##    b[2]     2.04   2.02 0.42 0.43   1.36   2.75 1.00      927      910
##    b[3]    -2.02  -2.01 0.52 0.52  -2.88  -1.18 1.00      762      945
##    s        1.02   1.00 0.15 0.14   0.80   1.29 1.00     2269     1757
##    y1[1]    2.18   2.20 1.09 1.07   0.37   3.95 1.00     1988     1906
##    y1[2]   -1.84  -1.85 1.11 1.07  -3.73  -0.03 1.00     1606     1894
##    y1[3]    0.17   0.19 1.09 1.08  -1.57   1.94 1.00     2122     2152
##    y1[4]    2.18   2.20 1.09 1.07   0.38   3.94 1.00     1944     1959
##    y1[5]    0.15   0.12 1.05 1.01  -1.57   1.90 1.00     1741     1980
##    y1[6]    2.23   2.22 1.07 1.03   0.52   3.99 1.00     1796     1852
##    y1[7]    0.17   0.16 1.07 1.02  -1.64   1.93 1.00     2066     1894
##    y1[8]    0.16   0.13 1.08 1.02  -1.58   1.96 1.00     1931     1847
##    y1[9]    0.13   0.11 1.06 1.04  -1.58   1.91 1.00     1719     1803
##    y1[10]   0.17   0.18 1.06 1.05  -1.60   1.84 1.00     1883     1715
##    y1[11]   0.13   0.11 1.08 1.08  -1.62   1.90 1.00     1943     1877
##    y1[12]   0.18   0.17 1.08 1.05  -1.53   2.03 1.00     1989     1811
##    y1[13]   0.19   0.18 1.08 1.09  -1.56   1.90 1.00     1890     1888
##    y1[14]   0.19   0.21 1.05 1.02  -1.59   1.86 1.00     1808     1934
##    y1[15]   0.15   0.15 1.07 1.02  -1.65   1.90 1.00     1776     1933
##    y1[16]   2.19   2.18 1.07 1.07   0.44   3.96 1.00     1841     1945
##    y1[17]   2.18   2.20 1.10 1.09   0.37   3.90 1.00     1973     1732
##    y1[18]   2.18   2.13 1.09 1.03   0.40   3.98 1.00     1698     1820
##    y1[19]   0.13   0.15 1.07 1.03  -1.68   1.86 1.00     1950     1885
##    y1[20]   2.23   2.23 1.05 1.02   0.51   3.92 1.00     1962     1781
##    y1[21]  -1.89  -1.91 1.14 1.08  -3.78  -0.01 1.00     1840     1915
##    y1[22]   2.18   2.21 1.09 1.07   0.34   3.87 1.00     1932     1934
##    y1[23]  -1.89  -1.90 1.09 1.06  -3.70  -0.08 1.00     1715     1942
##    y1[24]  -1.83  -1.82 1.09 1.05  -3.65  -0.01 1.00     1918     1981
##    y1[25]   2.22   2.22 1.08 1.07   0.53   4.00 1.00     1725     1878
##    y1[26]  -1.85  -1.84 1.12 1.09  -3.72   0.01 1.00     1882     1676
##    y1[27]   2.20   2.20 1.03 1.02   0.55   3.93 1.00     1990     1894
##    y1[28]   2.23   2.23 1.06 1.05   0.49   3.94 1.00     1974     1853
##    y1[29]   2.14   2.14 1.10 1.06   0.36   4.02 1.00     2038     2040
##    y1[30]  -1.89  -1.89 1.14 1.11  -3.68  -0.04 1.00     1832     1786
##    m1[1]    2.20   2.19 0.31 0.31   1.70   2.71 1.00     1410     1404
##    m1[2]   -1.86  -1.85 0.44 0.45  -2.57  -1.13 1.00     1004      924
##    m1[3]    0.16   0.17 0.29 0.27  -0.32   0.65 1.00      980     1057
##    m1[4]    2.20   2.19 0.31 0.31   1.70   2.71 1.00     1410     1404
##    m1[5]    0.16   0.17 0.29 0.27  -0.32   0.65 1.00      980     1057
##    m1[6]    2.20   2.19 0.31 0.31   1.70   2.71 1.00     1410     1404
##    m1[7]    0.16   0.17 0.29 0.27  -0.32   0.65 1.00      980     1057
##    m1[8]    0.16   0.17 0.29 0.27  -0.32   0.65 1.00      980     1057
##    m1[9]    0.16   0.17 0.29 0.27  -0.32   0.65 1.00      980     1057
##    m1[10]   0.16   0.17 0.29 0.27  -0.32   0.65 1.00      980     1057
##    m1[11]   0.16   0.17 0.29 0.27  -0.32   0.65 1.00      980     1057
##    m1[12]   0.16   0.17 0.29 0.27  -0.32   0.65 1.00      980     1057
##    m1[13]   0.16   0.17 0.29 0.27  -0.32   0.65 1.00      980     1057
##    m1[14]   0.16   0.17 0.29 0.27  -0.32   0.65 1.00      980     1057
##    m1[15]   0.16   0.17 0.29 0.27  -0.32   0.65 1.00      980     1057
##    m1[16]   2.20   2.19 0.31 0.31   1.70   2.71 1.00     1410     1404
##    m1[17]   2.20   2.19 0.31 0.31   1.70   2.71 1.00     1410     1404
##    m1[18]   2.20   2.19 0.31 0.31   1.70   2.71 1.00     1410     1404
##    m1[19]   0.16   0.17 0.29 0.27  -0.32   0.65 1.00      980     1057
##    m1[20]   2.20   2.19 0.31 0.31   1.70   2.71 1.00     1410     1404
##    m1[21]  -1.86  -1.85 0.44 0.45  -2.57  -1.13 1.00     1004      924
##    m1[22]   2.20   2.19 0.31 0.31   1.70   2.71 1.00     1410     1404
##    m1[23]  -1.86  -1.85 0.44 0.45  -2.57  -1.13 1.00     1004      924
##    m1[24]  -1.86  -1.85 0.44 0.45  -2.57  -1.13 1.00     1004      924
##    m1[25]   2.20   2.19 0.31 0.31   1.70   2.71 1.00     1410     1404
##    m1[26]  -1.86  -1.85 0.44 0.45  -2.57  -1.13 1.00     1004      924
##    m1[27]   2.20   2.19 0.31 0.31   1.70   2.71 1.00     1410     1404
##    m1[28]   2.20   2.19 0.31 0.31   1.70   2.71 1.00     1410     1404
##    m1[29]   2.20   2.19 0.31 0.31   1.70   2.71 1.00     1410     1404
##    m1[30]  -1.86  -1.85 0.44 0.45  -2.57  -1.13 1.00     1004      924

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

lv=c(0,1,2)
par(mfrow=c(1,1))
plot(tb$y,tb$y1,pch=lv[factor(tb$c)],xlab='obs.',ylab='prd.')

qplot(data=tb,y,y1,shape=c,size=I(2),xlab='obs.',ylab='prd.')

plot(tb$y,tb$y1,col=1+lv[factor(tb$c)],xlab='obs.',ylab='prd.')

qplot(data=tb,y,y1,col=c,xlab='obs.',ylab='prd.')

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(data=tb,1:n,y,col=c)+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

ANCOVA

tb=tibble(x=runif(n,0,9),c=sample(c('a','b','c'),n,replace=T),
          y=rnorm(n,x+(c=='b')*2-(c=='c')*2,1))

f=formula(y~x+c)
par(mfrow=c(1,1))
#plot(tb$x1,tb$y,pch=tb$c1)

lv=c(0,1,2)
plot(tb$x,tb$y,pch=lv[factor(tb$c)])

qplot(data=tb,x,y,shape=c,size=I(2))

plot(tb$x,tb$y,col=1+lv[factor(tb$c)])

qplot(data=tb,x,y,col=c)

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
## Initial log joint probability = -44312 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       42      -14.4242   0.000147739    0.00233842      0.6781      0.6781       94    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -14.42
##    b[1]      -0.04
##    b[2]       1.06
##    b[3]       1.38
##    b[4]      -2.43
##    s          0.98
##    y1[1]      5.22
##    y1[2]      0.54
##    y1[3]      0.88
##    y1[4]      4.00
##    y1[5]      4.53
##    y1[6]      1.56
##    y1[7]      9.10
##    y1[8]      4.33
##    y1[9]     -0.57
##    y1[10]     1.62
##    y1[11]     5.38
##    y1[12]     0.29
##    y1[13]     5.83
##    y1[14]     2.24
##    y1[15]     5.13
##    y1[16]    -0.73
##    y1[17]     2.41
##    y1[18]     8.63
##    y1[19]     6.02
##    y1[20]     4.08
##    y1[21]     3.69
##    y1[22]     7.15
##    y1[23]    -4.73
##    y1[24]     6.66
##    y1[25]     8.62
##    y1[26]     6.05
##    y1[27]    10.11
##    y1[28]     5.76
##    y1[29]     1.72
##    y1[30]     5.00
##    m1[1]      5.69
##    m1[2]      1.68
##    m1[3]      0.57
##    m1[4]      3.37
##    m1[5]      3.89
##    m1[6]     -0.14
##    m1[7]      9.54
##    m1[8]      4.88
##    m1[9]     -0.37
##    m1[10]     0.21
##    m1[11]     6.47
##    m1[12]     0.87
##    m1[13]     6.61
##    m1[14]     2.10
##    m1[15]     5.25
##    m1[16]     0.42
##    m1[17]     1.24
##    m1[18]     9.10
##    m1[19]     4.96
##    m1[20]     4.63
##    m1[21]     2.30
##    m1[22]     6.84
##    m1[23]    -2.07
##    m1[24]     7.42
##    m1[25]     7.68
##    m1[26]     7.63
##    m1[27]     8.21
##    m1[28]     4.64
##    m1[29]     2.54
##    m1[30]     4.95
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -17.35 -16.93 1.89 1.61 -20.93 -15.08 1.00      713      855
##    b[1]    -0.03  -0.03 0.50 0.47  -0.87   0.83 1.00      680      805
##    b[2]     1.05   1.05 0.09 0.09   0.91   1.19 1.00     1166     1111
##    b[3]     1.39   1.40 0.53 0.53   0.52   2.29 1.01      667      350
##    b[4]    -2.43  -2.45 0.48 0.46  -3.21  -1.65 1.00      714      919
##    s        1.12   1.10 0.17 0.16   0.88   1.42 1.00     1496     1120
##    y1[1]    5.68   5.69 1.22 1.23   3.68   7.64 1.00     1857     1864
##    y1[2]    1.73   1.72 1.17 1.15  -0.19   3.60 1.00     1689     1742
##    y1[3]    0.62   0.63 1.20 1.17  -1.34   2.54 1.00     1579     1847
##    y1[4]    3.42   3.42 1.24 1.20   1.34   5.47 1.00     1798     1768
##    y1[5]    3.84   3.85 1.17 1.15   1.89   5.72 1.00     1853     1856
##    y1[6]   -0.14  -0.15 1.18 1.15  -2.14   1.70 1.00     1884     1761
##    y1[7]    9.55   9.53 1.27 1.27   7.48  11.61 1.00     1819     1784
##    y1[8]    4.85   4.87 1.22 1.15   2.86   6.84 1.00     1809     1819
##    y1[9]   -0.38  -0.41 1.20 1.13  -2.33   1.61 1.00     1817     2050
##    y1[10]   0.20   0.19 1.14 1.13  -1.66   1.99 1.00     1612     1928
##    y1[11]   6.45   6.45 1.21 1.16   4.44   8.40 1.00     1957     2001
##    y1[12]   0.84   0.80 1.21 1.13  -1.13   2.89 1.00     2030     1775
##    y1[13]   6.62   6.64 1.18 1.08   4.64   8.52 1.00     1906     1931
##    y1[14]   2.07   2.07 1.19 1.16   0.14   4.04 1.00     1972     1904
##    y1[15]   5.20   5.21 1.20 1.12   3.18   7.12 1.00     1928     1859
##    y1[16]   0.41   0.40 1.21 1.19  -1.56   2.41 1.00     1826     1944
##    y1[17]   1.23   1.25 1.21 1.20  -0.76   3.21 1.00     1647     1933
##    y1[18]   9.11   9.13 1.21 1.19   7.13  11.08 1.00     1945     1839
##    y1[19]   4.98   5.00 1.17 1.10   3.01   6.86 1.00     2071     1807
##    y1[20]   4.61   4.59 1.18 1.15   2.73   6.57 1.00     1857     1708
##    y1[21]   2.31   2.31 1.19 1.18   0.39   4.20 1.00     1715     1704
##    y1[22]   6.84   6.85 1.20 1.15   4.99   8.83 1.00     1862     1967
##    y1[23]  -2.05  -2.07 1.23 1.20  -4.03   0.01 1.00     1855     1674
##    y1[24]   7.42   7.44 1.24 1.21   5.43   9.38 1.00     1814     1806
##    y1[25]   7.68   7.67 1.21 1.13   5.71   9.72 1.00     1602     1944
##    y1[26]   7.63   7.61 1.21 1.18   5.60   9.56 1.00     1936     1816
##    y1[27]   8.26   8.26 1.24 1.20   6.15  10.27 1.00     1328     1258
##    y1[28]   4.64   4.67 1.18 1.16   2.72   6.53 1.00     1768     1848
##    y1[29]   2.50   2.50 1.26 1.20   0.49   4.56 1.00     1743     1533
##    y1[30]   4.93   4.95 1.25 1.19   2.88   6.93 1.00     2054     1981
##    m1[1]    5.68   5.69 0.48 0.45   4.91   6.46 1.00     1180     1539
##    m1[2]    1.69   1.69 0.41 0.38   1.00   2.37 1.00      644      770
##    m1[3]    0.58   0.59 0.47 0.43  -0.21   1.37 1.00      664      829
##    m1[4]    3.39   3.39 0.50 0.49   2.56   4.21 1.00     1029      832
##    m1[5]    3.89   3.90 0.34 0.32   3.31   4.43 1.00      688      713
##    m1[6]   -0.14  -0.13 0.39 0.36  -0.79   0.51 1.00     1236     1344
##    m1[7]    9.54   9.53 0.51 0.48   8.70  10.39 1.00     1098     1159
##    m1[8]    4.86   4.88 0.44 0.42   4.15   5.58 1.00     1193     1552
##    m1[9]   -0.37  -0.37 0.40 0.37  -1.04   0.28 1.00     1242     1343
##    m1[10]   0.22   0.22 0.38 0.35  -0.42   0.84 1.00     1228     1308
##    m1[11]   6.48   6.48 0.44 0.42   5.73   7.19 1.00     1036      821
##    m1[12]   0.87   0.87 0.36 0.34   0.26   1.46 1.00     1221     1326
##    m1[13]   6.60   6.61 0.38 0.36   5.98   7.22 1.00     1023     1096
##    m1[14]   2.09   2.10 0.36 0.34   1.50   2.67 1.00     1245     1257
##    m1[15]   5.23   5.25 0.46 0.43   4.50   5.98 1.00     1191     1531
##    m1[16]   0.42   0.42 0.37 0.34  -0.20   1.03 1.00     1218     1364
##    m1[17]   1.25   1.25 0.43 0.40   0.52   1.97 1.00      650      823
##    m1[18]   9.08   9.08 0.50 0.48   8.28   9.92 1.00     1342     1137
##    m1[19]   4.95   4.96 0.34 0.32   4.36   5.49 1.00      783      856
##    m1[20]   4.63   4.63 0.34 0.32   4.05   5.17 1.00      748      782
##    m1[21]   2.31   2.31 0.38 0.36   1.67   2.94 1.00      641      711
##    m1[22]   6.83   6.84 0.39 0.38   6.20   7.46 1.00     1061     1188
##    m1[23]  -2.06  -2.07 0.47 0.45  -2.83  -1.30 1.00     1213     1175
##    m1[24]   7.41   7.41 0.41 0.41   6.75   8.08 1.00     1148     1318
##    m1[25]   7.69   7.68 0.45 0.43   6.93   8.43 1.00     1056      848
##    m1[26]   7.63   7.63 0.45 0.43   6.88   8.37 1.00     1055      839
##    m1[27]   8.22   8.21 0.46 0.44   7.44   8.98 1.00     1070      864
##    m1[28]   4.63   4.64 0.34 0.32   4.05   5.17 1.00      749      782
##    m1[29]   2.56   2.57 0.53 0.52   1.68   3.46 1.00     1039      940
##    m1[30]   4.93   4.95 0.44 0.42   4.22   5.65 1.00     1192     1532

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

lv=c(0,1,2)
par(mfrow=c(1,1))
plot(tb$y,tb$y1,pch=lv[factor(tb$c)],xlab='obs.',ylab='prd.')
abline(0,1)

plot(tb$y,tb$y1,col=1+lv[factor(tb$c)],xlab='obs.',ylab='prd.')
abline(0,1)

qplot(data=tb,y,y1,col=c,xlab='obs.',ylab='prd.')+
  geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(data=tb,1:n,y,col=c)+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

generalized linear regression

poisson regression

objective variable [0,Infinity), integer. variance of error is near to mean  
(normal linear regression, correlation is near to 1,-1, variance of error is constant)  

# of samples is N,  
log li=sum(bj*xji),j=[0,K],i=[1,N]  
yi~Po(li)  
 or  
li=sum(bj xji),j=[0,k]  
yi~Po(exp li)  

when xj -> xj +1, y -> y* exp bj   

ex6-1.stan

poisson regression

data{
  int N;
  int K;
  array[N] int y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
}
model{
  vector[N] l=X*b;
  y~poisson_log(l);
}
generated quantities{
  array[N] int y1;
  vector[N] l1=X*b;
  for(i in 1:N){
    y1[i]=poisson_log_rng(l1[i]);
  }
}
n=30
tb=tibble(x=runif(n,-1,1),c=sample(c('a','b'),n,replace=T),
          y=rpois(n,exp(x+(c=='b')*0.5)))
f=formula(y~x+c)
qplot(data=tb,x,y,col=c)

glm(f,tb,family='poisson')
## 
## Call:  glm(formula = f, family = "poisson", data = tb)
## 
## Coefficients:
## (Intercept)            x           cb  
##       0.035        0.822        0.344  
## 
## Degrees of Freedom: 29 Total (i.e. Null);  27 Residual
## Null Deviance:       36.7 
## Residual Deviance: 27.3  AIC: 94.7
X=model.matrix(f,tb)

data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mdl=cmdstan_model('./ex6-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -77.8209 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       12       -19.272   0.000160908   0.000572024      0.8673      0.8673       13    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -19.27
##    b[1]       0.03
##    b[2]       0.82
##    b[3]       0.34
##    y1[1]      2.00
##    y1[2]      1.00
##    y1[3]      2.00
##    y1[4]      2.00
##    y1[5]      5.00
##    y1[6]      0.00
##    y1[7]      0.00
##    y1[8]      0.00
##    y1[9]      3.00
##    y1[10]     1.00
##    y1[11]     0.00
##    y1[12]     1.00
##    y1[13]     1.00
##    y1[14]     2.00
##    y1[15]     1.00
##    y1[16]     0.00
##    y1[17]     2.00
##    y1[18]     0.00
##    y1[19]     2.00
##    y1[20]     0.00
##    y1[21]     2.00
##    y1[22]     1.00
##    y1[23]     3.00
##    y1[24]     0.00
##    y1[25]     1.00
##    y1[26]     5.00
##    y1[27]     2.00
##    y1[28]     2.00
##    y1[29]     3.00
##    y1[30]     4.00
##    l1[1]      0.04
##    l1[2]      0.37
##    l1[3]      0.86
##    l1[4]      0.90
##    l1[5]      0.71
##    l1[6]      0.70
##    l1[7]      0.57
##    l1[8]      0.59
##    l1[9]      0.78
##    l1[10]    -0.65
##    l1[11]    -0.59
##    l1[12]    -0.77
##    l1[13]     0.59
##    l1[14]     0.50
##    l1[15]     0.46
##    l1[16]     0.58
##    l1[17]    -0.08
##    l1[18]     0.15
##    l1[19]     0.82
##    l1[20]     0.52
##    l1[21]    -0.36
##    l1[22]     0.50
##    l1[23]     1.07
##    l1[24]    -0.20
##    l1[25]     0.26
##    l1[26]     0.65
##    l1[27]     0.80
##    l1[28]     0.70
##    l1[29]     1.12
##    l1[30]     1.08
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='l1',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -20.81 -20.49 1.29 1.03 -23.42 -19.45 1.01      816     1016
##    b[1]    -0.01   0.01 0.29 0.28  -0.54   0.43 1.02      669      826
##    b[2]     0.85   0.85 0.32 0.33   0.34   1.40 1.01      785      977
##    b[3]     0.36   0.34 0.30 0.29  -0.11   0.86 1.01      796      912
##    y1[1]    1.01   1.00 1.06 1.48   0.00   3.00 1.00     1892     1923
##    y1[2]    1.42   1.00 1.22 1.48   0.00   4.00 1.00     2005     1958
##    y1[3]    2.36   2.00 1.69 1.48   0.00   6.00 1.00     1796     1678
##    y1[4]    2.47   2.00 1.70 1.48   0.00   6.00 1.00     2219     1913
##    y1[5]    2.03   2.00 1.51 1.48   0.00   5.00 1.00     2057     1951
##    y1[6]    2.03   2.00 1.50 1.48   0.00   5.00 1.00     1772     1824
##    y1[7]    1.75   2.00 1.33 1.48   0.00   4.00 1.00     1961     1728
##    y1[8]    1.84   2.00 1.39 1.48   0.00   4.00 1.00     1884     1770
##    y1[9]    2.12   2.00 1.53 1.48   0.00   5.00 1.00     1868     2035
##    y1[10]   0.54   0.00 0.78 0.00   0.00   2.00 1.00     1753     1847
##    y1[11]   0.56   0.00 0.78 0.00   0.00   2.00 1.00     1917     1745
##    y1[12]   0.47   0.00 0.72 0.00   0.00   2.00 1.00     1660     1694
##    y1[13]   1.81   2.00 1.41 1.48   0.00   4.00 1.00     1680     1700
##    y1[14]   1.66   1.00 1.34 1.48   0.00   4.00 1.00     1705     1801
##    y1[15]   1.58   1.00 1.28 1.48   0.00   4.00 1.00     2003     1845
##    y1[16]   1.72   2.00 1.32 1.48   0.00   4.00 1.00     1904     1904
##    y1[17]   0.95   1.00 1.03 1.48   0.00   3.00 1.00     1860     1703
##    y1[18]   1.17   1.00 1.17 1.48   0.00   3.00 1.00     1982     2021
##    y1[19]   2.28   2.00 1.61 1.48   0.00   5.00 1.00     1749     1564
##    y1[20]   1.70   1.00 1.37 1.48   0.00   4.00 1.00     1806     1915
##    y1[21]   0.69   0.00 0.88 0.00   0.00   2.00 1.00     1852     1902
##    y1[22]   1.61   1.00 1.34 1.48   0.00   4.00 1.00     1743     1931
##    y1[23]   2.98   3.00 1.85 1.48   0.00   6.00 1.00     1670     1960
##    y1[24]   0.85   1.00 0.99 1.48   0.00   3.00 1.00     1782     1936
##    y1[25]   1.28   1.00 1.14 1.48   0.00   3.00 1.00     1842     1796
##    y1[26]   1.96   2.00 1.48 1.48   0.00   5.00 1.00     2084     2058
##    y1[27]   2.28   2.00 1.64 1.48   0.00   5.00 1.00     1770     2010
##    y1[28]   2.00   2.00 1.50 1.48   0.00   5.00 1.00     1790     1856
##    y1[29]   3.10   3.00 1.88 1.48   0.00   7.00 1.00     1674     1492
##    y1[30]   2.96   3.00 1.84 1.48   0.00   6.00 1.00     1933     1985
##    l1[1]   -0.01   0.00 0.31 0.31  -0.54   0.48 1.01     1053     1124
##    l1[2]    0.34   0.34 0.22 0.22  -0.03   0.69 1.00     1450     1474
##    l1[3]    0.84   0.84 0.19 0.18   0.51   1.15 1.00     1468     1381
##    l1[4]    0.89   0.89 0.19 0.18   0.55   1.21 1.00     1367     1407
##    l1[5]    0.68   0.70 0.25 0.25   0.28   1.07 1.00     1203     1153
##    l1[6]    0.68   0.68 0.18 0.17   0.39   0.97 1.00     1789     1387
##    l1[7]    0.54   0.54 0.19 0.18   0.24   0.84 1.00     1781     1652
##    l1[8]    0.57   0.58 0.24 0.23   0.18   0.93 1.01     1088      998
##    l1[9]    0.76   0.76 0.18 0.17   0.45   1.06 1.00     1655     1408
##    l1[10]  -0.72  -0.69 0.50 0.49  -1.60   0.06 1.02      645      884
##    l1[11]  -0.66  -0.63 0.48 0.47  -1.51   0.08 1.02      642      877
##    l1[12]  -0.84  -0.82 0.54 0.53  -1.80   0.00 1.02      651      863
##    l1[13]   0.57   0.58 0.24 0.23   0.18   0.93 1.01     1087     1009
##    l1[14]   0.47   0.48 0.23 0.23   0.07   0.83 1.01      969      926
##    l1[15]   0.42   0.43 0.20 0.20   0.08   0.75 1.00     1589     1434
##    l1[16]   0.55   0.55 0.18 0.18   0.25   0.86 1.00     1792     1544
##    l1[17]  -0.13  -0.12 0.35 0.36  -0.72   0.42 1.01      998      983
##    l1[18]   0.11   0.12 0.28 0.28  -0.37   0.55 1.01     1150     1147
##    l1[19]   0.80   0.80 0.18 0.17   0.48   1.11 1.00     1558     1387
##    l1[20]   0.49   0.50 0.23 0.23   0.09   0.85 1.01      994      947
##    l1[21]  -0.42  -0.40 0.41 0.40  -1.13   0.20 1.02      638      936
##    l1[22]   0.47   0.49 0.23 0.23   0.08   0.84 1.01      976      945
##    l1[23]   1.06   1.06 0.23 0.22   0.67   1.42 1.00     1099     1194
##    l1[24]  -0.25  -0.23 0.39 0.40  -0.91   0.36 1.01      960     1032
##    l1[25]   0.22   0.23 0.25 0.25  -0.20   0.61 1.00     1285     1327
##    l1[26]   0.63   0.65 0.24 0.24   0.23   1.00 1.00     1153     1024
##    l1[27]   0.78   0.79 0.26 0.26   0.34   1.19 1.00     1235     1160
##    l1[28]   0.68   0.69 0.25 0.25   0.27   1.06 1.00     1196     1153
##    l1[29]   1.11   1.12 0.24 0.24   0.71   1.51 1.00     1042     1102
##    l1[30]   1.07   1.07 0.23 0.23   0.68   1.44 1.00     1101     1193

y1=mcmc$draws('y1')
smy=summary(y1)

table(tb$y,smy$median)
##    
##     0 1 2 3
##   0 2 2 2 0
##   1 1 3 4 0
##   2 1 5 3 1
##   3 0 0 3 0
##   4 0 0 0 1
##   5 0 0 1 1
cat('\n')
table(tb$y,smy$median,tb$c)
## , ,  = a
## 
##    
##     0 1 2 3
##   0 2 0 1 0
##   1 1 0 2 0
##   2 1 3 1 0
##   3 0 0 2 0
##   4 0 0 0 0
##   5 0 0 0 0
## 
## , ,  = b
## 
##    
##     0 1 2 3
##   0 0 2 1 0
##   1 0 3 2 0
##   2 0 2 2 1
##   3 0 0 1 0
##   4 0 0 0 1
##   5 0 0 1 1
par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')

map=c()
for(i in 1:n){
  a=table(y1[,,i])
  map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)
##    map
##     0 1 2 3
##   0 3 3 0 0
##   1 4 3 1 0
##   2 1 7 1 1
##   3 0 1 2 0
##   4 0 0 0 1
##   5 0 0 2 0
cat('\n')
table(tb$y,map,tb$c)
## , ,  = a
## 
##    map
##     0 1 2 3
##   0 2 1 0 0
##   1 1 2 0 0
##   2 1 4 0 0
##   3 0 1 1 0
##   4 0 0 0 0
##   5 0 0 0 0
## 
## , ,  = b
## 
##    map
##     0 1 2 3
##   0 1 2 0 0
##   1 3 1 1 0
##   2 0 3 1 1
##   3 0 0 1 0
##   4 0 0 0 1
##   5 0 0 2 0
par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')

logistic regression

# of samples is N,  
objective variable 0/1 binary  
  
probability of incident pi[0,1]  
logit pi=log(pi/ 1-pi)=sum(bj*xji),j=[0,K],i=[1,N] (-Infinity, Infinity)  

yi~Ber(pi), 0/1 binary  

odds(x)=p(x)/ 1-p(x), probablity of incident / probablity of no incident  
odds ratio(x0,x1)=odds(x1)/odd(x0)  
  
when xj -> xj +1, odds ratio -> odds ratio *exp bj  

ex6-2.stan

logistic regression

data{
  int N;
  int K;
  array[N] int y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
}
model{
  vector[N] z=X*b;
  y~bernoulli_logit(z);
}
generated quantities{
  array[N] int y1;
  vector[N] z1=X*b;
  for(i in 1:N){
    y1[i]=bernoulli_rng(inv_logit(z1[i]));
  }
}
n=30
x=runif(n,-1,1)
c=sample(c('a','b'),n,replace=T)
z=x+(c=='b')
y=rbinom(n,1,1/(1+exp(-z)))
tb=tibble(x=x,c=c,y=y)

f=formula(y~x+c)
qplot(data=tb,x,y,col=c)

glm(f,tb,family='binomial') # it can caluculte when all trials are once
## 
## Call:  glm(formula = f, family = "binomial", data = tb)
## 
## Coefficients:
## (Intercept)            x           cb  
##      -0.890        0.645        0.937  
## 
## Degrees of Freedom: 29 Total (i.e. Null);  27 Residual
## Null Deviance:       41.1 
## Residual Deviance: 38.6  AIC: 44.6
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mdl=cmdstan_model('./ex6-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -37.8535 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       12      -19.2851   0.000136794   3.65965e-05      0.9577      0.9577       15    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -19.29
##    b[1]      -0.89
##    b[2]       0.64
##    b[3]       0.94
##    y1[1]      0.00
##    y1[2]      0.00
##    y1[3]      1.00
##    y1[4]      0.00
##    y1[5]      1.00
##    y1[6]      0.00
##    y1[7]      0.00
##    y1[8]      1.00
##    y1[9]      1.00
##    y1[10]     1.00
##    y1[11]     1.00
##    y1[12]     1.00
##    y1[13]     0.00
##    y1[14]     0.00
##    y1[15]     1.00
##    y1[16]     0.00
##    y1[17]     1.00
##    y1[18]     0.00
##    y1[19]     1.00
##    y1[20]     0.00
##    y1[21]     1.00
##    y1[22]     0.00
##    y1[23]     0.00
##    y1[24]     0.00
##    y1[25]     1.00
##    y1[26]     1.00
##    y1[27]     1.00
##    y1[28]     0.00
##    y1[29]     1.00
##    y1[30]     0.00
##    z1[1]     -0.94
##    z1[2]      0.54
##    z1[3]     -1.18
##    z1[4]     -0.83
##    z1[5]      0.45
##    z1[6]     -0.33
##    z1[7]     -0.26
##    z1[8]      0.49
##    z1[9]     -0.17
##    z1[10]     0.13
##    z1[11]    -0.81
##    z1[12]    -0.31
##    z1[13]     0.56
##    z1[14]     0.35
##    z1[15]     0.01
##    z1[16]    -1.43
##    z1[17]     0.29
##    z1[18]    -0.03
##    z1[19]     0.01
##    z1[20]     0.67
##    z1[21]     0.07
##    z1[22]    -0.57
##    z1[23]     0.13
##    z1[24]    -1.20
##    z1[25]    -1.11
##    z1[26]    -0.59
##    z1[27]    -0.45
##    z1[28]    -0.56
##    z1[29]    -1.48
##    z1[30]    -0.37
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='z1',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -20.86 -20.50 1.30 1.04 -23.43 -19.46 1.00      762      991
##    b[1]    -1.00  -0.96 0.67 0.66  -2.17   0.03 1.00      957      881
##    b[2]     0.75   0.73 0.71 0.68  -0.38   1.96 1.00     1249     1041
##    b[3]     1.03   1.00 0.82 0.80  -0.25   2.41 1.01      900      737
##    y1[1]    0.26   0.00 0.44 0.00   0.00   1.00 1.00     1794       NA
##    y1[2]    0.63   1.00 0.48 0.00   0.00   1.00 1.00     1912       NA
##    y1[3]    0.22   0.00 0.42 0.00   0.00   1.00 1.00     1718       NA
##    y1[4]    0.30   0.00 0.46 0.00   0.00   1.00 1.00     1988       NA
##    y1[5]    0.63   1.00 0.48 0.00   0.00   1.00 1.00     1672       NA
##    y1[6]    0.43   0.00 0.50 0.00   0.00   1.00 1.00     2053       NA
##    y1[7]    0.44   0.00 0.50 0.00   0.00   1.00 1.00     1810       NA
##    y1[8]    0.61   1.00 0.49 0.00   0.00   1.00 1.00     2153       NA
##    y1[9]    0.43   0.00 0.50 0.00   0.00   1.00 1.00     1927       NA
##    y1[10]   0.55   1.00 0.50 0.00   0.00   1.00 1.00     2113       NA
##    y1[11]   0.31   0.00 0.46 0.00   0.00   1.00 1.00     2022       NA
##    y1[12]   0.40   0.00 0.49 0.00   0.00   1.00 1.00     1811       NA
##    y1[13]   0.63   1.00 0.48 0.00   0.00   1.00 1.00     2035       NA
##    y1[14]   0.60   1.00 0.49 0.00   0.00   1.00 1.00     1904       NA
##    y1[15]   0.51   1.00 0.50 0.00   0.00   1.00 1.00     1949       NA
##    y1[16]   0.21   0.00 0.41 0.00   0.00   1.00 1.00     1942       NA
##    y1[17]   0.57   1.00 0.49 0.00   0.00   1.00 1.00     1908       NA
##    y1[18]   0.48   0.00 0.50 0.00   0.00   1.00 1.00     1845       NA
##    y1[19]   0.49   0.00 0.50 0.00   0.00   1.00 1.00     1938       NA
##    y1[20]   0.64   1.00 0.48 0.00   0.00   1.00 1.00     2129       NA
##    y1[21]   0.51   1.00 0.50 0.00   0.00   1.00 1.00     1741       NA
##    y1[22]   0.34   0.00 0.47 0.00   0.00   1.00 1.00     1913       NA
##    y1[23]   0.53   1.00 0.50 0.00   0.00   1.00 1.00     2129       NA
##    y1[24]   0.23   0.00 0.42 0.00   0.00   1.00 1.00     1893       NA
##    y1[25]   0.25   0.00 0.43 0.00   0.00   1.00 1.00     1771       NA
##    y1[26]   0.35   0.00 0.48 0.00   0.00   1.00 1.00     1735       NA
##    y1[27]   0.40   0.00 0.49 0.00   0.00   1.00 1.00     1574       NA
##    y1[28]   0.34   0.00 0.47 0.00   0.00   1.00 1.00     1938       NA
##    y1[29]   0.20   0.00 0.40 0.00   0.00   1.00 1.00     1877       NA
##    y1[30]   0.40   0.00 0.49 0.00   0.00   1.00 1.00     1903       NA
##    z1[1]   -1.07  -1.03 0.69 0.67  -2.30  -0.02 1.00      949      881
##    z1[2]    0.59   0.56 0.69 0.66  -0.48   1.78 1.00     1591     1499
##    z1[3]   -1.35  -1.29 0.81 0.79  -2.77  -0.14 1.00      957      919
##    z1[4]   -0.93  -0.89 0.66 0.64  -2.06   0.08 1.00      971      934
##    z1[5]    0.49   0.47 0.63 0.60  -0.49   1.58 1.00     1687     1404
##    z1[6]   -0.35  -0.34 0.79 0.79  -1.68   0.93 1.00     1233     1218
##    z1[7]   -0.27  -0.25 0.84 0.82  -1.69   1.11 1.00     1269     1219
##    z1[8]    0.54   0.52 0.66 0.63  -0.47   1.68 1.00     1635     1414
##    z1[9]   -0.23  -0.22 0.61 0.59  -1.27   0.78 1.00     1830     1392
##    z1[10]   0.12   0.11 0.52 0.51  -0.76   0.98 1.00     2093     1166
##    z1[11]  -0.91  -0.87 0.66 0.64  -2.05   0.10 1.00      975      916
##    z1[12]  -0.39  -0.38 0.70 0.68  -1.57   0.77 1.00     1677     1310
##    z1[13]   0.63   0.59 0.71 0.68  -0.48   1.86 1.00     1567     1573
##    z1[14]   0.38   0.36 0.58 0.54  -0.55   1.39 1.00     1823     1269
##    z1[15]  -0.02  -0.02 0.53 0.53  -0.92   0.87 1.00     2049     1127
##    z1[16]  -1.63  -1.58 0.99 0.96  -3.34  -0.16 1.00      990      891
##    z1[17]   0.30   0.29 0.55 0.53  -0.59   1.25 1.00     1923     1176
##    z1[18]  -0.07  -0.06 0.55 0.54  -1.00   0.84 1.00     2013     1263
##    z1[19]  -0.02  -0.02 0.53 0.53  -0.92   0.87 1.00     2052     1127
##    z1[20]   0.75   0.72 0.79 0.76  -0.50   2.11 1.00     1490     1683
##    z1[21]   0.05   0.04 0.52 0.52  -0.86   0.91 1.00     2094     1266
##    z1[22]  -0.63  -0.60 0.68 0.66  -1.81   0.42 1.00     1093     1086
##    z1[23]   0.12   0.12 0.52 0.51  -0.76   0.99 1.00     2091     1164
##    z1[24]  -1.36  -1.31 0.82 0.79  -2.82  -0.13 1.00      958      925
##    z1[25]  -1.26  -1.21 0.77 0.75  -2.61  -0.10 1.00      950      919
##    z1[26]  -0.72  -0.70 0.94 0.89  -2.29   0.81 1.00     1474     1285
##    z1[27]  -0.49  -0.46 0.73 0.71  -1.75   0.65 1.00     1168     1054
##    z1[28]  -0.69  -0.67 0.91 0.87  -2.21   0.81 1.00     1488     1286
##    z1[29]  -1.69  -1.63 1.04 1.00  -3.48  -0.13 1.00      998      891
##    z1[30]  -0.39  -0.37 0.77 0.77  -1.70   0.84 1.00     1217     1143

y1=mcmc$draws('y1')
smy=summary(y1)

table(tb$y,smy$median)
##    
##      0  1
##   0 13  4
##   1  6  7
cat('\n')
table(tb$y,smy$median,tb$c)
## , ,  = a
## 
##    
##     0 1
##   0 9 0
##   1 4 0
## 
## , ,  = b
## 
##    
##     0 1
##   0 4 4
##   1 2 7
par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')

map=c()
for(i in 1:n){
  a=table(y1[,,i])
  map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)
##    map
##      0  1
##   0 13  4
##   1  6  7
cat('\n')
table(tb$y,map,tb$c)
## , ,  = a
## 
##    map
##     0 1
##   0 9 0
##   1 4 0
## 
## , ,  = b
## 
##    map
##     0 1
##   0 4 4
##   1 2 7
par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')

multi logistic regression

# of samples is N,  
# of trials of a sample i is mi,  
objective variable [0,n], integer  
  
probability of incident pi[0,1]  
logit pi=log(pi/ 1-pi)=sum(bj*xji),j=[0,K],i=[1,N] (-Infinity, Infinity)  

yi~B(mi, pi), # of acutual incident  

odds(x)=p(x)/ 1-p(x), probablity of incident / probablity of no incident  
odds ratio(x0,x1)=odds(x1)/odd(x0)  
  
when xj -> xj +1, odds ratio -> odds ratio *exp bj  

ex6-3.stan

multi logistic regression

data{
  int N;
  int K;
  array[N] int m;
  array[N] int y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
}
model{
  vector[N] z=X*b;
  y~binomial_logit(m,z);
}
generated quantities{
  array[N] int y1;
  vector[N] z1=X*b;
  for(i in 1:N){
    y1[i]=binomial_rng(m[i],inv_logit(z1[i]));
  }
}
n=30
m=floor(runif(n,1,10)) # trials are varying (1,10)
x=runif(n,-1,1)
c=sample(c('a','b'),n,replace=T)
z=x+(c=='b')
y=rbinom(n,m,1/(1+exp(-z)))
tb=tibble(x=x,c=c,y=y)

f=formula(y~x+c)
qplot(data=tb,x,y,col=c)

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X,m=m)

mdl=cmdstan_model('./ex6-3.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -148.389 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        8      -95.2055    0.00130982    0.00101506           1           1       10    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -95.21
##    b[1]      -0.23
##    b[2]       0.83
##    b[3]       1.20
##    y1[1]      1.00
##    y1[2]      2.00
##    y1[3]      2.00
##    y1[4]      3.00
##    y1[5]      1.00
##    y1[6]      7.00
##    y1[7]      2.00
##    y1[8]      5.00
##    y1[9]      1.00
##    y1[10]     1.00
##    y1[11]     4.00
##    y1[12]     2.00
##    y1[13]     3.00
##    y1[14]     6.00
##    y1[15]     2.00
##    y1[16]     3.00
##    y1[17]     1.00
##    y1[18]     4.00
##    y1[19]     3.00
##    y1[20]     1.00
##    y1[21]     0.00
##    y1[22]     5.00
##    y1[23]     5.00
##    y1[24]     3.00
##    y1[25]     1.00
##    y1[26]     5.00
##    y1[27]     3.00
##    y1[28]     6.00
##    y1[29]     1.00
##    y1[30]     5.00
##    z1[1]      0.43
##    z1[2]      1.57
##    z1[3]     -0.95
##    z1[4]      1.01
##    z1[5]      0.20
##    z1[6]      0.99
##    z1[7]      0.35
##    z1[8]     -0.14
##    z1[9]      1.03
##    z1[10]     1.66
##    z1[11]     0.04
##    z1[12]     0.35
##    z1[13]    -0.92
##    z1[14]     1.19
##    z1[15]     0.34
##    z1[16]    -0.23
##    z1[17]     1.57
##    z1[18]    -0.38
##    z1[19]    -1.00
##    z1[20]    -0.44
##    z1[21]    -0.86
##    z1[22]     0.51
##    z1[23]     1.49
##    z1[24]    -0.99
##    z1[25]    -0.63
##    z1[26]     1.22
##    z1[27]    -0.07
##    z1[28]     0.41
##    z1[29]     0.53
##    z1[30]     0.58
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='z1',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -96.69 -96.38 1.23 0.97 -99.18 -95.39 1.00      891     1092
##    b[1]    -0.25  -0.24 0.25 0.26  -0.68   0.16 1.00     1447     1449
##    b[2]     0.84   0.84 0.35 0.35   0.29   1.42 1.00     1147     1135
##    b[3]     1.24   1.25 0.33 0.34   0.70   1.79 1.00     1148     1229
##    y1[1]    1.80   2.00 0.87 1.48   0.00   3.00 1.00     1932       NA
##    y1[2]    1.65   2.00 0.55 0.00   1.00   2.00 1.00     2012       NA
##    y1[3]    2.51   2.00 1.47 1.48   0.00   5.00 1.00     1843     1976
##    y1[4]    2.95   3.00 0.89 1.48   1.00   4.00 1.00     1989       NA
##    y1[5]    1.11   1.00 0.70 1.48   0.00   2.00 1.00     2093       NA
##    y1[6]    6.58   7.00 1.42 1.48   4.00   9.00 1.00     2008       NA
##    y1[7]    3.53   4.00 1.29 1.48   1.00   6.00 1.00     1948       NA
##    y1[8]    2.74   3.00 1.26 1.48   1.00   5.00 1.00     1986     2026
##    y1[9]    3.66   4.00 1.02 1.48   2.00   5.00 1.00     1824       NA
##    y1[10]   0.84   1.00 0.37 0.00   0.00   1.00 1.00     2113       NA
##    y1[11]   4.03   4.00 1.55 1.48   2.00   7.00 1.00     1786     1956
##    y1[12]   1.16   1.00 0.72 1.48   0.00   2.00 1.00     1924       NA
##    y1[13]   1.38   1.00 1.04 1.48   0.00   3.00 1.00     2178     2014
##    y1[14]   6.87   7.00 1.34 1.48   5.00   9.00 1.00     1854       NA
##    y1[15]   2.89   3.00 1.18 1.48   1.00   5.00 1.00     1754       NA
##    y1[16]   3.54   4.00 1.46 1.48   1.00   6.00 1.00     1912     2029
##    y1[17]   0.83   1.00 0.38 0.00   0.00   1.00 1.00     2110       NA
##    y1[18]   1.97   2.00 1.15 1.48   0.00   4.00 1.00     1882     1900
##    y1[19]   1.87   2.00 1.23 1.48   0.00   4.00 1.00     1740     1763
##    y1[20]   1.59   2.00 0.98 1.48   0.00   3.00 1.00     2039     1671
##    y1[21]   0.30   0.00 0.46 0.00   0.00   1.00 1.00     1736       NA
##    y1[22]   4.38   4.00 1.34 1.48   2.00   6.00 1.00     1919     2058
##    y1[23]   7.33   8.00 1.23 1.48   5.00   9.00 1.00     1810       NA
##    y1[24]   2.20   2.00 1.33 1.48   0.00   4.00 1.00     1717     1801
##    y1[25]   0.66   1.00 0.68 1.48   0.00   2.00 1.00     2209       NA
##    y1[26]   3.89   4.00 0.97 1.48   2.00   5.00 1.00     1897       NA
##    y1[27]   1.93   2.00 1.02 1.48   0.00   4.00 1.00     1944       NA
##    y1[28]   4.88   5.00 1.50 1.48   2.00   7.00 1.00     1861     1940
##    y1[29]   0.65   1.00 0.48 0.00   0.00   1.00 1.00     1992       NA
##    y1[30]   5.13   5.00 1.40 1.48   3.00   7.00 1.00     2012     1939
##    z1[1]    0.44   0.45 0.29 0.28  -0.04   0.91 1.00     1702     1534
##    z1[2]    1.60   1.58 0.40 0.40   0.95   2.26 1.00     1230     1327
##    z1[3]   -0.98  -0.96 0.33 0.33  -1.53  -0.45 1.00     1231     1410
##    z1[4]    1.03   1.03 0.25 0.25   0.61   1.45 1.00     1658     1653
##    z1[5]    0.22   0.22 0.35 0.34  -0.35   0.79 1.00     1565     1569
##    z1[6]    1.02   1.02 0.25 0.25   0.60   1.44 1.00     1674     1692
##    z1[7]    0.37   0.37 0.31 0.31  -0.13   0.87 1.00     1689     1622
##    z1[8]   -0.16  -0.16 0.27 0.27  -0.61   0.27 1.00     1423     1448
##    z1[9]    1.05   1.05 0.26 0.25   0.63   1.48 1.00     1630     1772
##    z1[10]   1.69   1.68 0.43 0.43   1.00   2.41 1.00     1206     1307
##    z1[11]   0.02   0.03 0.31 0.31  -0.50   0.54 1.00     1371     1354
##    z1[12]   0.34   0.34 0.40 0.40  -0.32   1.00 1.00     1290     1358
##    z1[13]  -0.95  -0.93 0.32 0.32  -1.49  -0.43 1.00     1238     1401
##    z1[14]   1.21   1.21 0.29 0.28   0.75   1.70 1.00     1442     1470
##    z1[15]   0.33   0.33 0.40 0.40  -0.33   0.98 1.00     1292     1358
##    z1[16]  -0.25  -0.24 0.25 0.26  -0.67   0.17 1.00     1448     1449
##    z1[17]   1.60   1.59 0.40 0.40   0.95   2.26 1.00     1229     1327
##    z1[18]  -0.40  -0.40 0.24 0.24  -0.80  -0.01 1.00     1448     1546
##    z1[19]  -1.03  -1.01 0.35 0.35  -1.60  -0.48 1.00     1219     1378
##    z1[20]  -0.46  -0.46 0.24 0.24  -0.86  -0.07 1.00     1436     1475
##    z1[21]  -0.88  -0.87 0.31 0.31  -1.40  -0.40 1.00     1259     1444
##    z1[22]   0.53   0.52 0.27 0.27   0.08   0.96 1.00     1768     1626
##    z1[23]   1.52   1.51 0.38 0.37   0.92   2.15 1.00     1255     1424
##    z1[24]  -1.02  -1.01 0.34 0.35  -1.59  -0.47 1.00     1220     1378
##    z1[25]  -0.66  -0.65 0.26 0.26  -1.10  -0.25 1.00     1356     1467
##    z1[26]   1.25   1.24 0.30 0.29   0.77   1.74 1.00     1413     1447
##    z1[27]  -0.08  -0.08 0.28 0.29  -0.57   0.37 1.00     1402     1536
##    z1[28]   0.43   0.43 0.29 0.29  -0.05   0.91 1.00     1692     1585
##    z1[29]   0.55   0.55 0.27 0.26   0.11   0.98 1.00     1787     1742
##    z1[30]   0.60   0.60 0.26 0.26   0.17   1.01 1.00     1822     1759

y1=mcmc$draws('y1')
smy=summary(y1)

table(tb$y,smy$median)
##    
##     0 1 2 3 4 5 7 8
##   0 1 2 0 0 0 0 0 0
##   1 0 4 3 0 0 0 0 0
##   2 0 0 4 0 0 0 0 0
##   3 0 1 1 3 1 1 0 0
##   4 0 0 0 0 1 0 0 0
##   5 0 0 0 0 4 1 1 0
##   6 0 0 0 0 0 0 1 0
##   9 0 0 0 0 0 0 0 1
cat('\n')
table(tb$y,smy$median,tb$c)
## , ,  = a
## 
##    
##     0 1 2 3 4 5 7 8
##   0 1 1 0 0 0 0 0 0
##   1 0 1 3 0 0 0 0 0
##   2 0 0 3 0 0 0 0 0
##   3 0 1 0 2 0 0 0 0
##   4 0 0 0 0 0 0 0 0
##   5 0 0 0 0 2 0 0 0
##   6 0 0 0 0 0 0 0 0
##   9 0 0 0 0 0 0 0 0
## 
## , ,  = b
## 
##    
##     0 1 2 3 4 5 7 8
##   0 0 1 0 0 0 0 0 0
##   1 0 3 0 0 0 0 0 0
##   2 0 0 1 0 0 0 0 0
##   3 0 0 1 1 1 1 0 0
##   4 0 0 0 0 1 0 0 0
##   5 0 0 0 0 2 1 1 0
##   6 0 0 0 0 0 0 1 0
##   9 0 0 0 0 0 0 0 1
par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')

map=c()
for(i in 1:n){
  a=table(y1[,,i])
  map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)
##    map
##     0 1 2 3 4 5 7 8
##   0 1 2 0 0 0 0 0 0
##   1 1 4 2 0 0 0 0 0
##   2 0 0 4 0 0 0 0 0
##   3 0 1 1 3 1 1 0 0
##   4 0 0 0 0 0 1 0 0
##   5 0 0 0 0 4 1 1 0
##   6 0 0 0 0 0 0 1 0
##   9 0 0 0 0 0 0 0 1
cat('\n')
table(tb$y,map,tb$c)
## , ,  = a
## 
##    map
##     0 1 2 3 4 5 7 8
##   0 1 1 0 0 0 0 0 0
##   1 1 1 2 0 0 0 0 0
##   2 0 0 3 0 0 0 0 0
##   3 0 1 0 2 0 0 0 0
##   4 0 0 0 0 0 0 0 0
##   5 0 0 0 0 2 0 0 0
##   6 0 0 0 0 0 0 0 0
##   9 0 0 0 0 0 0 0 0
## 
## , ,  = b
## 
##    map
##     0 1 2 3 4 5 7 8
##   0 0 1 0 0 0 0 0 0
##   1 0 3 0 0 0 0 0 0
##   2 0 0 1 0 0 0 0 0
##   3 0 0 1 1 1 1 0 0
##   4 0 0 0 0 0 1 0 0
##   5 0 0 0 0 2 1 1 0
##   6 0 0 0 0 0 0 1 0
##   9 0 0 0 0 0 0 0 1
par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')

interaction of variable

n=50
mdl=cmdstan_model('./ex5.stan')

tb=tibble(x=runif(n,-3,3),
          ca=sample(c('a','b'),n,replace=T),cb=sample(c('a','b'),n,replace=T),
          y=rnorm(n,x+(ca=='b')+(cb=='b')-(ca=='b')*(cb=='b'),1))

grid.arrange(qplot(data=tb,x,y,col=ca),
             qplot(data=tb,x,y,col=cb),ncol=2)

fn=function(f){
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
mle

mcmc=goMCMC(mdl,data)

seeMCMC(mcmc,exc='m',ch=F)

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

grid.arrange(
  qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
grid.arrange(
  qplot(data=tb,1:n,y,col=ca)+
    geom_point(aes(x=1:n,y=y1),tb,col='black')+
    geom_line(aes(x=1:n,y=l5),tb,col='gray')+
    geom_line(aes(x=1:n,y=u5),tb,col='gray'),
  qplot(data=tb,1:n,y,col=cb)+
    geom_point(aes(x=1:n,y=y1),tb,col='black')+
    geom_line(aes(x=1:n,y=l5),tb,col='gray')+
    geom_line(aes(x=1:n,y=u5),tb,col='gray'),
  ncol=2
)
}


f0=formula(y~x+ca)
f1=formula(y~x+ca+cb)
f2=formula(y~x+ca*cb)

fn(f0)
## Initial log joint probability = -966.17 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       12      -25.9661   0.000395456    0.00220611           1           1       17    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -27.95 -27.62 1.45 1.21 -30.84 -26.32 1.00      812     1070
##    b[1]     0.35   0.36 0.21 0.21   0.01   0.67 1.00      598      538
##    b[2]     0.88   0.87 0.09 0.09   0.72   1.02 1.00     2570     1543
##    b[3]     0.33   0.32 0.29 0.28  -0.12   0.82 1.00      566      601
##    s        1.08   1.07 0.12 0.11   0.91   1.29 1.00     2142     1403
##    y1[1]    0.70   0.70 1.12 1.09  -1.15   2.50 1.00     2036     1943
##    y1[2]    1.08   1.08 1.11 1.08  -0.77   2.92 1.00     1968     2026
##    y1[3]   -2.03  -2.01 1.11 1.12  -3.84  -0.22 1.00     1620     1751
##    y1[4]    1.25   1.24 1.12 1.08  -0.58   3.06 1.00     1934     2050
##    y1[5]    2.57   2.57 1.13 1.11   0.72   4.50 1.00     1818     1881
##    y1[6]    1.97   1.96 1.09 1.07   0.17   3.83 1.00     2038     1900
##    y1[7]    0.28   0.30 1.11 1.10  -1.50   2.12 1.00     1838     1855
##    y1[8]    2.32   2.31 1.12 1.14   0.53   4.19 1.00     1966     1933
##    y1[9]   -1.40  -1.40 1.14 1.12  -3.28   0.46 1.00     1814     1914
##    y1[10]  -2.28  -2.30 1.16 1.15  -4.18  -0.35 1.00     1818     1959
##    y1[11]   1.58   1.60 1.14 1.15  -0.32   3.44 1.00     1896     1882
##    y1[12]   0.13   0.14 1.11 1.07  -1.69   1.96 1.00     1990     1952
##    y1[13]  -0.89  -0.89 1.09 1.05  -2.67   0.93 1.00     1870     1899
##    y1[14]   0.23   0.22 1.11 1.09  -1.56   2.08 1.00     2098     1933
##    y1[15]   2.24   2.26 1.14 1.16   0.41   4.06 1.00     1913     2103
##    y1[16]  -0.03  -0.03 1.11 1.11  -1.83   1.88 1.00     1843     1822
##    y1[17]  -0.83  -0.81 1.09 1.05  -2.64   0.92 1.00     2094     2032
##    y1[18]   0.23   0.23 1.11 1.09  -1.54   2.07 1.00     1981     1892
##    y1[19]   0.56   0.56 1.11 1.11  -1.24   2.46 1.00     1918     1773
##    y1[20]   1.69   1.69 1.10 1.09  -0.15   3.55 1.00     1740     1884
##    y1[21]  -0.29  -0.28 1.08 1.06  -2.08   1.48 1.00     1947     1857
##    y1[22]   2.97   2.98 1.15 1.11   1.10   4.87 1.00     1897     1822
##    y1[23]  -0.58  -0.57 1.11 1.12  -2.42   1.23 1.00     1916     2016
##    y1[24]   2.99   3.01 1.16 1.12   1.08   4.93 1.00     2089     1818
##    y1[25]  -1.38  -1.36 1.13 1.09  -3.29   0.42 1.00     1940     1958
##    y1[26]   1.84   1.85 1.12 1.10   0.00   3.67 1.00     1965     1703
##    y1[27]   1.62   1.62 1.12 1.09  -0.26   3.43 1.00     1843     1893
##    y1[28]  -1.64  -1.62 1.12 1.08  -3.49   0.19 1.00     1826     1916
##    y1[29]   2.93   2.93 1.13 1.08   1.10   4.77 1.00     1985     1786
##    y1[30]  -1.65  -1.62 1.10 1.10  -3.54   0.09 1.00     1772     1840
##    y1[31]   1.31   1.28 1.10 1.07  -0.49   3.11 1.00     1869     1861
##    y1[32]   0.35   0.37 1.10 1.09  -1.50   2.18 1.00     1834     1961
##    y1[33]  -0.75  -0.76 1.11 1.08  -2.60   1.06 1.00     2027     1910
##    y1[34]  -0.74  -0.75 1.09 1.05  -2.52   1.07 1.00     1976     1972
##    y1[35]  -1.53  -1.53 1.10 1.10  -3.29   0.28 1.00     1925     1863
##    y1[36]  -0.08  -0.10 1.09 1.07  -1.92   1.73 1.00     1936     2002
##    y1[37]  -0.42  -0.40 1.10 1.09  -2.26   1.33 1.00     1888     1879
##    y1[38]   1.11   1.14 1.09 1.07  -0.64   2.92 1.00     1920     1891
##    y1[39]  -0.55  -0.56 1.11 1.08  -2.35   1.28 1.00     2117     1848
##    y1[40]  -1.12  -1.11 1.10 1.05  -2.97   0.69 1.00     1688     1934
##    y1[41]  -2.25  -2.27 1.12 1.13  -4.06  -0.44 1.00     1784     1821
##    y1[42]   0.31   0.31 1.12 1.07  -1.54   2.18 1.00     1921     1658
##    y1[43]   1.62   1.60 1.09 1.04  -0.15   3.44 1.00     1891     1866
##    y1[44]   2.48   2.47 1.14 1.16   0.65   4.40 1.00     1845     1521
##    y1[45]   1.36   1.39 1.11 1.12  -0.49   3.08 1.00     2037     1893
##    y1[46]  -1.48  -1.48 1.14 1.11  -3.32   0.37 1.00     2132     1813
##    y1[47]  -1.13  -1.14 1.11 1.08  -2.93   0.72 1.00     1973     1969
##    y1[48]  -1.96  -1.98 1.14 1.13  -3.80  -0.14 1.00     1782     1893
##    y1[49]  -0.61  -0.60 1.09 1.09  -2.43   1.17 1.00     1826     1965
##    y1[50]   0.78   0.81 1.10 1.05  -1.09   2.57 1.00     1735     1776
##    m1[1]    0.73   0.73 0.21 0.21   0.39   1.07 1.00     1259     1333
##    m1[2]    1.05   1.05 0.22 0.21   0.69   1.40 1.00     1300     1286
##    m1[3]   -1.98  -1.98 0.29 0.28  -2.46  -1.49 1.00     1358     1357
##    m1[4]    1.26   1.27 0.24 0.24   0.86   1.65 1.00      631      759
##    m1[5]    2.60   2.59 0.34 0.33   2.03   3.15 1.00      849     1120
##    m1[6]    2.01   2.01 0.26 0.26   1.59   2.45 1.00     1535     1224
##    m1[7]    0.28   0.28 0.21 0.21  -0.08   0.63 1.00     1262     1279
##    m1[8]    2.33   2.33 0.28 0.28   1.87   2.80 1.00     1620     1220
##    m1[9]   -1.36  -1.36 0.29 0.29  -1.84  -0.90 1.00     1752     1445
##    m1[10]  -2.24  -2.24 0.31 0.30  -2.75  -1.73 1.00     1515     1581
##    m1[11]   1.56   1.57 0.26 0.26   1.13   1.98 1.00      665      978
##    m1[12]   0.13   0.14 0.21 0.21  -0.20   0.45 1.00      611      546
##    m1[13]  -0.92  -0.91 0.26 0.26  -1.35  -0.49 1.00     1511     1446
##    m1[14]   0.23   0.23 0.21 0.21  -0.13   0.57 1.00     1263     1279
##    m1[15]   2.28   2.28 0.32 0.30   1.76   2.80 1.00      787      990
##    m1[16]  -0.01  -0.01 0.22 0.21  -0.37   0.35 1.00     1286     1436
##    m1[17]  -0.86  -0.87 0.22 0.22  -1.23  -0.49 1.00      822      955
##    m1[18]   0.19   0.20 0.21 0.21  -0.14   0.51 1.00      607      502
##    m1[19]   0.55   0.56 0.21 0.21   0.20   0.88 1.00      595      627
##    m1[20]   1.69   1.69 0.27 0.26   1.24   2.13 1.00      685      957
##    m1[21]  -0.29  -0.29 0.23 0.22  -0.67   0.09 1.00     1305     1368
##    m1[22]   2.99   2.99 0.33 0.32   2.44   3.55 1.00     1787     1184
##    m1[23]  -0.59  -0.59 0.21 0.21  -0.94  -0.24 1.00      737      728
##    m1[24]   2.99   2.99 0.33 0.32   2.44   3.54 1.00     1786     1184
##    m1[25]  -1.40  -1.41 0.25 0.25  -1.82  -0.98 1.00     1052     1282
##    m1[26]   1.82   1.83 0.25 0.25   1.41   2.25 1.00     1492     1247
##    m1[27]   1.60   1.61 0.26 0.26   1.16   2.03 1.00      671      894
##    m1[28]  -1.64  -1.64 0.31 0.32  -2.14  -1.13 1.00     1835     1435
##    m1[29]   2.92   2.92 0.33 0.32   2.38   3.47 1.00     1774     1198
##    m1[30]  -1.69  -1.70 0.27 0.26  -2.14  -1.24 1.00     1204     1282
##    m1[31]   1.31   1.31 0.25 0.24   0.90   1.70 1.00      635      815
##    m1[32]   0.34   0.34 0.21 0.21  -0.02   0.68 1.00     1255     1222
##    m1[33]  -0.74  -0.74 0.25 0.25  -1.15  -0.33 1.00     1431     1403
##    m1[34]  -0.72  -0.72 0.22 0.22  -1.07  -0.36 1.00      776      844
##    m1[35]  -1.49  -1.48 0.30 0.30  -1.97  -1.00 1.00     1791     1413
##    m1[36]  -0.08  -0.07 0.20 0.21  -0.41   0.24 1.00      635      567
##    m1[37]  -0.41  -0.41 0.23 0.23  -0.80  -0.03 1.00     1331     1369
##    m1[38]   1.07   1.07 0.22 0.21   0.71   1.42 1.00     1303     1286
##    m1[39]  -0.57  -0.57 0.24 0.24  -0.96  -0.18 1.00     1369     1330
##    m1[40]  -1.13  -1.14 0.24 0.23  -1.53  -0.74 1.00      929     1075
##    m1[41]  -2.25  -2.25 0.31 0.30  -2.76  -1.74 1.00     1520     1556
##    m1[42]   0.32   0.32 0.21 0.21  -0.04   0.66 1.00     1256     1222
##    m1[43]   1.60   1.61 0.24 0.23   1.21   2.00 1.00     1423     1196
##    m1[44]   2.47   2.47 0.33 0.31   1.92   3.01 1.00      824     1081
##    m1[45]   1.40   1.40 0.23 0.22   1.02   1.77 1.00     1368     1347
##    m1[46]  -1.44  -1.45 0.25 0.25  -1.87  -1.01 1.00     1069     1282
##    m1[47]  -1.11  -1.10 0.27 0.27  -1.56  -0.67 1.00     1617     1499
##    m1[48]  -1.95  -1.95 0.29 0.27  -2.43  -1.47 1.00     1341     1320
##    m1[49]  -0.62  -0.62 0.21 0.21  -0.97  -0.27 1.00      747      728
##    m1[50]   0.81   0.82 0.22 0.22   0.45   1.15 1.00      601      665

fn(f1)
## Initial log joint probability = -60.4972 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       18      -25.2382   0.000130537     0.0018028           1           1       22    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -27.78 -27.47 1.62 1.51 -30.75 -25.82 1.00      840     1101
##    b[1]     0.13   0.14 0.27 0.28  -0.31   0.56 1.01      708     1084
##    b[2]     0.88   0.88 0.09 0.09   0.74   1.03 1.00     2780     1882
##    b[3]     0.46   0.45 0.32 0.31  -0.06   0.97 1.01      693     1001
##    b[4]     0.37   0.38 0.32 0.32  -0.16   0.90 1.01      820      916
##    s        1.07   1.07 0.12 0.11   0.90   1.28 1.00     2431     1653
##    y1[1]    0.63   0.63 1.09 1.05  -1.22   2.40 1.00     2009     1931
##    y1[2]    1.38   1.38 1.14 1.11  -0.49   3.24 1.00     1937     1971
##    y1[3]   -1.83  -1.81 1.12 1.08  -3.75  -0.02 1.00     1979     1729
##    y1[4]    1.44   1.47 1.16 1.12  -0.50   3.30 1.00     1923     1750
##    y1[5]    2.79   2.81 1.11 1.11   1.01   4.59 1.00     1696     2023
##    y1[6]    1.92   1.90 1.12 1.06   0.06   3.80 1.00     2029     1854
##    y1[7]    0.21   0.20 1.11 1.10  -1.64   2.04 1.00     1961     1969
##    y1[8]    2.29   2.28 1.15 1.15   0.42   4.18 1.00     1786     1748
##    y1[9]   -1.45  -1.44 1.14 1.13  -3.36   0.45 1.00     1875     1880
##    y1[10]  -2.10  -2.11 1.12 1.08  -4.01  -0.34 1.00     2036     1877
##    y1[11]   1.37   1.37 1.13 1.07  -0.53   3.25 1.00     1698     1927
##    y1[12]   0.29   0.27 1.08 1.07  -1.45   2.10 1.00     2120     2011
##    y1[13]  -1.03  -1.04 1.10 1.09  -2.86   0.75 1.00     1730     1745
##    y1[14]   0.50   0.49 1.12 1.08  -1.32   2.39 1.00     1901     1955
##    y1[15]   2.07   2.04 1.17 1.15   0.15   3.92 1.00     1567     1807
##    y1[16]  -0.10  -0.09 1.10 1.09  -1.92   1.71 1.00     1819     1966
##    y1[17]  -1.10  -1.10 1.12 1.13  -2.92   0.74 1.00     1860     1858
##    y1[18]   0.37   0.32 1.11 1.11  -1.44   2.23 1.00     1527     1776
##    y1[19]   0.71   0.70 1.09 1.10  -1.07   2.49 1.00     1908     1878
##    y1[20]   1.46   1.46 1.11 1.08  -0.38   3.28 1.00     1842     1767
##    y1[21]  -0.05  -0.07 1.14 1.08  -1.95   1.84 1.00     1854     1693
##    y1[22]   2.96   2.95 1.15 1.11   1.12   4.79 1.00     1719     1971
##    y1[23]  -0.43  -0.45 1.12 1.09  -2.25   1.43 1.00     1829     1933
##    y1[24]   2.90   2.90 1.13 1.12   1.01   4.77 1.00     1927     1952
##    y1[25]  -1.25  -1.21 1.14 1.08  -3.16   0.59 1.00     1882     1790
##    y1[26]   2.12   2.13 1.16 1.13   0.20   4.00 1.00     1921     1929
##    y1[27]   1.40   1.41 1.11 1.15  -0.44   3.19 1.00     1742     1881
##    y1[28]  -1.77  -1.74 1.11 1.11  -3.57   0.06 1.00     1672     1863
##    y1[29]   2.87   2.85 1.16 1.17   0.95   4.76 1.00     2035     1787
##    y1[30]  -1.92  -1.95 1.13 1.13  -3.75  -0.05 1.00     1723     1955
##    y1[31]   1.48   1.48 1.12 1.07  -0.35   3.33 1.00     1830     1831
##    y1[32]   0.21   0.22 1.12 1.13  -1.59   2.09 1.00     2025     1936
##    y1[33]  -0.87  -0.86 1.10 1.09  -2.65   0.94 1.00     2032     1976
##    y1[34]  -0.62  -0.61 1.08 1.07  -2.41   1.14 1.00     2168     2004
##    y1[35]  -1.61  -1.62 1.14 1.15  -3.44   0.41 1.00     2018     2035
##    y1[36]  -0.31  -0.30 1.10 1.08  -2.14   1.49 1.00     2014     2015
##    y1[37]  -0.14  -0.16 1.15 1.10  -2.06   1.78 1.00     1937     1821
##    y1[38]   0.96   0.97 1.13 1.10  -0.92   2.82 1.00     1963     1895
##    y1[39]  -0.65  -0.67 1.11 1.06  -2.45   1.24 1.00     1757     1881
##    y1[40]  -1.42  -1.41 1.11 1.09  -3.26   0.47 1.00     1742     1712
##    y1[41]  -2.10  -2.10 1.12 1.10  -3.93  -0.25 1.00     1932     1899
##    y1[42]   0.58   0.56 1.12 1.08  -1.24   2.47 1.00     1995     1713
##    y1[43]   1.50   1.48 1.11 1.14  -0.34   3.34 1.00     1978     2013
##    y1[44]   2.64   2.64 1.10 1.11   0.81   4.43 1.00     1975     1822
##    y1[45]   1.31   1.33 1.10 1.09  -0.47   3.06 1.00     1936     1893
##    y1[46]  -1.30  -1.29 1.09 1.09  -3.10   0.41 1.01     1582     1574
##    y1[47]  -1.23  -1.25 1.12 1.13  -3.04   0.62 1.00     2106     1903
##    y1[48]  -2.17  -2.15 1.17 1.13  -4.07  -0.31 1.00     1921     1972
##    y1[49]  -0.46  -0.45 1.12 1.10  -2.36   1.41 1.00     2031     1867
##    y1[50]   0.62   0.64 1.12 1.10  -1.20   2.47 1.00     1918     2019
##    m1[1]    0.64   0.64 0.23 0.23   0.26   1.04 1.00     1681     1603
##    m1[2]    1.34   1.35 0.33 0.32   0.78   1.87 1.01      894     1005
##    m1[3]   -1.84  -1.84 0.32 0.32  -2.37  -1.32 1.00     1896     1802
##    m1[4]    1.43   1.42 0.29 0.30   0.94   1.90 1.00     1484     1569
##    m1[5]    2.77   2.76 0.38 0.39   2.16   3.39 1.00     1754     1712
##    m1[6]    1.93   1.93 0.28 0.27   1.49   2.39 1.00     2168     1535
##    m1[7]    0.19   0.19 0.23 0.23  -0.19   0.58 1.00     1596     1606
##    m1[8]    2.26   2.26 0.29 0.29   1.78   2.74 1.00     2295     1753
##    m1[9]   -1.47  -1.47 0.31 0.30  -1.98  -0.98 1.00     1694     1749
##    m1[10]  -2.11  -2.11 0.34 0.33  -2.67  -1.56 1.00     1988     1852
##    m1[11]   1.35   1.35 0.31 0.30   0.85   1.84 1.00      868     1310
##    m1[12]   0.29   0.29 0.26 0.25  -0.14   0.71 1.00     1398     1294
##    m1[13]  -1.02  -1.02 0.28 0.28  -1.48  -0.55 1.00     1629     1780
##    m1[14]   0.51   0.52 0.32 0.31  -0.03   1.03 1.01      873     1062
##    m1[15]   2.08   2.08 0.35 0.35   1.51   2.65 1.00     1048     1542
##    m1[16]  -0.11  -0.10 0.24 0.24  -0.50   0.29 1.00     1572     1537
##    m1[17]  -1.09  -1.09 0.29 0.29  -1.55  -0.61 1.01      760     1062
##    m1[18]   0.35   0.35 0.26 0.25  -0.09   0.77 1.00     1398     1293
##    m1[19]   0.71   0.71 0.27 0.26   0.26   1.14 1.00     1406     1299
##    m1[20]   1.48   1.48 0.31 0.31   0.97   1.98 1.00      897     1367
##    m1[21]  -0.01   0.00 0.32 0.32  -0.56   0.51 1.01      912     1132
##    m1[22]   2.92   2.93 0.34 0.35   2.37   3.48 1.00     2503     1819
##    m1[23]  -0.44  -0.43 0.26 0.26  -0.86  -0.02 1.00     1483     1530
##    m1[24]   2.92   2.92 0.34 0.35   2.36   3.47 1.00     2502     1819
##    m1[25]  -1.26  -1.25 0.29 0.28  -1.73  -0.79 1.00     1700     1750
##    m1[26]   2.12   2.12 0.36 0.34   1.53   2.70 1.01      994     1146
##    m1[27]   1.39   1.39 0.31 0.31   0.88   1.88 1.00      877     1329
##    m1[28]  -1.75  -1.74 0.33 0.32  -2.29  -1.23 1.00     1741     1814
##    m1[29]   2.85   2.86 0.34 0.34   2.31   3.40 1.00     2484     1820
##    m1[30]  -1.92  -1.92 0.33 0.33  -2.46  -1.38 1.01      902     1126
##    m1[31]   1.47   1.47 0.30 0.30   0.98   1.95 1.00     1493     1568
##    m1[32]   0.24   0.24 0.23 0.23  -0.14   0.63 1.00     1604     1606
##    m1[33]  -0.84  -0.84 0.27 0.26  -1.29  -0.39 1.00     1615     1703
##    m1[34]  -0.57  -0.57 0.27 0.26  -1.00  -0.14 1.00     1509     1601
##    m1[35]  -1.59  -1.59 0.31 0.31  -2.11  -1.09 1.00     1713     1816
##    m1[36]  -0.30  -0.29 0.27 0.27  -0.74   0.14 1.01      702     1073
##    m1[37]  -0.14  -0.12 0.33 0.32  -0.69   0.39 1.01      925     1217
##    m1[38]   0.98   0.98 0.24 0.23   0.60   1.39 1.00     1796     1603
##    m1[39]  -0.67  -0.67 0.26 0.26  -1.11  -0.24 1.00     1593     1637
##    m1[40]  -1.36  -1.36 0.30 0.30  -1.84  -0.85 1.01      803     1118
##    m1[41]  -2.11  -2.11 0.34 0.33  -2.67  -1.57 1.00     1991     1852
##    m1[42]   0.60   0.61 0.32 0.31   0.06   1.13 1.01      869     1072
##    m1[43]   1.52   1.53 0.26 0.24   1.10   1.97 1.00     1992     1651
##    m1[44]   2.64   2.63 0.37 0.38   2.04   3.24 1.00     1728     1713
##    m1[45]   1.31   1.32 0.25 0.24   0.91   1.74 1.00     1932     1536
##    m1[46]  -1.30  -1.29 0.29 0.28  -1.78  -0.83 1.00     1712     1719
##    m1[47]  -1.21  -1.21 0.29 0.28  -1.70  -0.73 1.00     1661     1749
##    m1[48]  -2.18  -2.18 0.35 0.35  -2.75  -1.61 1.01      960     1211
##    m1[49]  -0.47  -0.47 0.26 0.26  -0.90  -0.05 1.00     1488     1530
##    m1[50]   0.60   0.60 0.28 0.28   0.16   1.03 1.00      743     1166

fn(f2)
## Initial log joint probability = -285.114 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       37      -23.3564   6.67346e-05   0.000830854      0.3772      0.9909       46    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -26.61 -26.29 1.81 1.72 -30.07 -24.24 1.00      687     1057
##    b[1]    -0.18  -0.18 0.33 0.31  -0.71   0.36 1.00      635      996
##    b[2]     0.89   0.89 0.09 0.09   0.75   1.03 1.00     2164     1312
##    b[3]     0.94   0.94 0.41 0.41   0.27   1.61 1.00      529      796
##    b[4]     0.89   0.90 0.43 0.40   0.20   1.62 1.00      561      935
##    b[5]    -1.21  -1.21 0.64 0.66  -2.21  -0.16 1.01      414      563
##    s        1.05   1.04 0.12 0.12   0.88   1.26 1.00     2322     1499
##    y1[1]    0.85   0.85 1.09 1.07  -0.96   2.65 1.00     1936     1813
##    y1[2]    0.81   0.82 1.16 1.12  -1.08   2.74 1.00     1647     1819
##    y1[3]   -1.65  -1.68 1.08 1.10  -3.46   0.08 1.00     1970     2042
##    y1[4]    1.64   1.63 1.11 1.08  -0.18   3.48 1.00     2018     1833
##    y1[5]    3.02   3.02 1.09 1.06   1.26   4.84 1.00     1900     1860
##    y1[6]    2.13   2.14 1.10 1.10   0.32   3.90 1.00     2106     1542
##    y1[7]    0.41   0.41 1.07 1.04  -1.28   2.16 1.00     1960     2014
##    y1[8]    2.43   2.45 1.08 1.05   0.64   4.20 1.00     1983     1930
##    y1[9]   -1.33  -1.36 1.07 1.06  -3.12   0.47 1.00     1869     1932
##    y1[10]  -1.95  -1.95 1.10 1.10  -3.70  -0.11 1.00     1791     1453
##    y1[11]   1.02   1.03 1.10 1.11  -0.75   2.88 1.00     1822     1733
##    y1[12]   0.49   0.49 1.09 1.10  -1.26   2.26 1.00     1731     1794
##    y1[13]  -0.89  -0.89 1.11 1.10  -2.72   0.92 1.00     1759     1621
##    y1[14]  -0.01  -0.01 1.13 1.11  -1.88   1.84 1.00     1663     1903
##    y1[15]   1.81   1.83 1.14 1.15  -0.11   3.69 1.00     2007     1847
##    y1[16]   0.04   0.04 1.06 1.07  -1.73   1.75 1.00     1945     1879
##    y1[17]  -1.46  -1.43 1.10 1.11  -3.33   0.30 1.00     1769     1793
##    y1[18]   0.55   0.52 1.09 1.11  -1.21   2.35 1.00     2039     1954
##    y1[19]   0.92   0.93 1.10 1.12  -0.84   2.67 1.00     1883     1973
##    y1[20]   1.21   1.19 1.12 1.12  -0.58   3.03 1.00     1627     1777
##    y1[21]  -0.56  -0.48 1.13 1.08  -2.41   1.24 1.00     1948     1750
##    y1[22]   3.14   3.11 1.11 1.08   1.36   4.93 1.00     1963     2010
##    y1[23]  -0.21  -0.23 1.07 1.06  -1.97   1.59 1.00     1750     1643
##    y1[24]   3.14   3.16 1.10 1.10   1.35   4.99 1.00     1915     1746
##    y1[25]  -1.05  -1.04 1.10 1.06  -2.92   0.74 1.00     1870     1892
##    y1[26]   1.66   1.68 1.14 1.10  -0.26   3.53 1.00     1791     2012
##    y1[27]   1.09   1.08 1.10 1.11  -0.72   2.88 1.00     1545     1708
##    y1[28]  -1.63  -1.65 1.10 1.08  -3.37   0.18 1.00     2063     1845
##    y1[29]   3.05   3.08 1.13 1.11   1.20   4.83 1.00     1945     1984
##    y1[30]  -2.27  -2.23 1.11 1.11  -4.11  -0.46 1.00     1642     1668
##    y1[31]   1.66   1.67 1.12 1.09  -0.14   3.48 1.00     1855     1738
##    y1[32]   0.42   0.40 1.07 1.10  -1.29   2.15 1.00     2092     2055
##    y1[33]  -0.70  -0.69 1.13 1.14  -2.53   1.11 1.00     2057     1918
##    y1[34]  -0.38  -0.39 1.09 1.05  -2.13   1.41 1.00     2115     1904
##    y1[35]  -1.44  -1.41 1.12 1.11  -3.28   0.37 1.00     1966     1867
##    y1[36]  -0.60  -0.63 1.11 1.09  -2.38   1.26 1.00     1606     1919
##    y1[37]  -0.66  -0.67 1.15 1.12  -2.48   1.27 1.00     2062     1899
##    y1[38]   1.19   1.19 1.09 1.06  -0.58   2.95 1.00     2183     2027
##    y1[39]  -0.52  -0.52 1.09 1.07  -2.29   1.28 1.00     2024     1855
##    y1[40]  -1.69  -1.68 1.13 1.14  -3.59   0.12 1.00     1528     1913
##    y1[41]  -1.91  -1.92 1.10 1.05  -3.71  -0.06 1.00     1979     1972
##    y1[42]   0.07   0.08 1.15 1.11  -1.85   2.01 1.00     1498     1631
##    y1[43]   1.77   1.77 1.07 1.07   0.01   3.53 1.00     1967     1974
##    y1[44]   2.92   2.91 1.12 1.08   1.12   4.70 1.00     1828     2013
##    y1[45]   1.46   1.44 1.07 1.05  -0.30   3.24 1.00     1860     1792
##    y1[46]  -1.12  -1.10 1.09 1.05  -2.90   0.70 1.00     1943     1739
##    y1[47]  -1.11  -1.11 1.10 1.08  -2.93   0.69 1.00     1996     1889
##    y1[48]  -2.53  -2.53 1.12 1.07  -4.41  -0.67 1.00     1584     1666
##    y1[49]  -0.31  -0.29 1.08 1.09  -2.10   1.45 1.00     1634     1719
##    y1[50]   0.30   0.31 1.10 1.11  -1.54   2.08 1.00     1834     2109
##    m1[1]    0.82   0.82 0.25 0.24   0.40   1.23 1.00     1488     1411
##    m1[2]    0.83   0.83 0.44 0.42   0.14   1.55 1.00     1115     1072
##    m1[3]   -1.65  -1.66 0.33 0.34  -2.16  -1.10 1.00     1645     1161
##    m1[4]    1.64   1.65 0.32 0.33   1.11   2.16 1.00     1355     1497
##    m1[5]    3.01   3.01 0.40 0.40   2.33   3.65 1.00     1514     1627
##    m1[6]    2.12   2.12 0.29 0.28   1.65   2.60 1.00     1572     1550
##    m1[7]    0.36   0.36 0.25 0.24  -0.05   0.77 1.00     1494     1402
##    m1[8]    2.45   2.45 0.31 0.30   1.94   2.95 1.00     1615     1452
##    m1[9]   -1.32  -1.32 0.31 0.31  -1.81  -0.79 1.00     1737     1441
##    m1[10]  -1.92  -1.93 0.34 0.35  -2.46  -1.35 1.00     1713     1280
##    m1[11]   1.05   1.06 0.35 0.34   0.47   1.64 1.00      754     1116
##    m1[12]   0.50   0.50 0.28 0.29   0.03   0.96 1.00     1308     1245
##    m1[13]  -0.86  -0.86 0.29 0.29  -1.31  -0.37 1.00     1647     1510
##    m1[14]  -0.01  -0.01 0.43 0.42  -0.71   0.73 1.00     1052     1158
##    m1[15]   1.79   1.80 0.39 0.37   1.17   2.46 1.00      885     1278
##    m1[16]   0.06   0.06 0.25 0.25  -0.36   0.48 1.00     1513     1423
##    m1[17]  -1.41  -1.42 0.35 0.33  -1.96  -0.85 1.00      651     1133
##    m1[18]   0.56   0.56 0.28 0.29   0.08   1.01 1.00     1309     1246
##    m1[19]   0.92   0.92 0.29 0.29   0.44   1.39 1.00     1312     1214
##    m1[20]   1.19   1.20 0.36 0.35   0.60   1.79 1.00      776     1116
##    m1[21]  -0.54  -0.54 0.44 0.42  -1.26   0.22 1.00     1028     1219
##    m1[22]   3.12   3.12 0.35 0.34   2.56   3.70 1.00     1707     1373
##    m1[23]  -0.24  -0.24 0.28 0.28  -0.69   0.23 1.00     1356     1302
##    m1[24]   3.12   3.12 0.35 0.34   2.56   3.69 1.00     1707     1373
##    m1[25]  -1.07  -1.07 0.30 0.31  -1.55  -0.56 1.00     1504     1253
##    m1[26]   1.62   1.61 0.45 0.45   0.89   2.36 1.00     1207     1194
##    m1[27]   1.09   1.10 0.36 0.34   0.51   1.68 1.00      760     1117
##    m1[28]  -1.60  -1.61 0.33 0.33  -2.12  -1.04 1.00     1791     1394
##    m1[29]   3.05   3.05 0.35 0.33   2.49   3.62 1.00     1698     1389
##    m1[30]  -2.25  -2.26 0.38 0.36  -2.85  -1.61 1.00      726     1170
##    m1[31]   1.69   1.70 0.32 0.33   1.15   2.21 1.00     1359     1490
##    m1[32]   0.41   0.41 0.25 0.24   0.00   0.82 1.00     1494     1362
##    m1[33]  -0.68  -0.68 0.28 0.28  -1.13  -0.21 1.00     1613     1487
##    m1[34]  -0.37  -0.37 0.28 0.28  -0.82   0.10 1.00     1377     1159
##    m1[35]  -1.44  -1.45 0.32 0.31  -1.95  -0.90 1.00     1766     1394
##    m1[36]  -0.61  -0.61 0.33 0.31  -1.14  -0.07 1.00      627     1001
##    m1[37]  -0.66  -0.67 0.44 0.42  -1.39   0.10 1.00     1026     1235
##    m1[38]   1.16   1.16 0.26 0.25   0.74   1.57 1.00     1489     1456
##    m1[39]  -0.51  -0.51 0.27 0.27  -0.94  -0.06 1.00     1581     1406
##    m1[40]  -1.69  -1.70 0.35 0.34  -2.25  -1.10 1.00      670     1135
##    m1[41]  -1.93  -1.94 0.34 0.35  -2.47  -1.36 1.00     1715     1280
##    m1[42]   0.08   0.08 0.43 0.42  -0.62   0.82 1.00     1057     1143
##    m1[43]   1.71   1.71 0.27 0.26   1.25   2.15 1.00     1513     1444
##    m1[44]   2.87   2.88 0.39 0.39   2.21   3.50 1.00     1497     1656
##    m1[45]   1.50   1.50 0.26 0.25   1.06   1.94 1.00     1499     1367
##    m1[46]  -1.11  -1.11 0.30 0.31  -1.59  -0.59 1.00     1514     1271
##    m1[47]  -1.06  -1.06 0.30 0.30  -1.52  -0.55 1.00     1687     1511
##    m1[48]  -2.52  -2.52 0.39 0.38  -3.14  -1.85 1.00      756     1150
##    m1[49]  -0.27  -0.28 0.28 0.29  -0.72   0.20 1.00     1361     1219
##    m1[50]   0.29   0.30 0.33 0.32  -0.26   0.83 1.00      665      988

tb=tibble(xa=runif(n,-2,2),xb=runif(n,-2,2),
          ca=sample(c('a','b'),n,replace=T),cb=sample(c('a','b'),n,replace=T),
          y=rnorm(n,xa+xb-xa*xb+(ca=='b')*2+(cb=='b')*2-(ca=='b')*(cb=='b'),1))

grid.arrange(qplot(data=tb,xa,y,col=ca),
             qplot(data=tb,xa,y,col=cb),
             qplot(data=tb,xb,y,col=ca),
             qplot(data=tb,xb,y,col=cb))

fn=function(f){
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
mle

mcmc=goMCMC(mdl,data)

seeMCMC(mcmc,exc='m',ch=F)

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

grid.arrange(
  qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
grid.arrange(
  qplot(data=tb,1:n,y,col=ca)+
    geom_point(aes(x=1:n,y=y1),tb,col='black')+
    geom_line(aes(x=1:n,y=l5),tb,col='gray')+
    geom_line(aes(x=1:n,y=u5),tb,col='gray'),
  qplot(data=tb,1:n,y,col=cb)+
    geom_point(aes(x=1:n,y=y1),tb,col='black')+
    geom_line(aes(x=1:n,y=l5),tb,col='gray')+
    geom_line(aes(x=1:n,y=u5),tb,col='gray'),
  ncol=2
)
}


f0=formula(y~xa+xb+ca+cb)
f1=formula(y~xa+xb+ca*cb)
f2=formula(y~xa*xb+ca*cb)

fn(f0)
## Initial log joint probability = -83.506 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       21      -49.0092   0.000177649    0.00121943           1           1       23    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -51.86 -51.46 1.96 1.71 -55.60 -49.40 1.00      774     1028
##    b[1]     1.03   1.03 0.41 0.39   0.35   1.69 1.01      866      902
##    b[2]     1.03   1.03 0.26 0.25   0.61   1.46 1.00     2134     1495
##    b[3]     1.01   1.00 0.25 0.24   0.60   1.41 1.00     2136     1270
##    b[4]     0.04   0.05 0.50 0.46  -0.78   0.86 1.00      929      939
##    b[5]     2.07   2.06 0.52 0.51   1.23   2.94 1.01      832      879
##    s        1.76   1.74 0.20 0.19   1.47   2.11 1.01     1913     1489
##    y1[1]    0.74   0.75 1.89 1.84  -2.36   3.87 1.00     1876     1971
##    y1[2]   -0.33  -0.29 1.84 1.86  -3.33   2.60 1.00     1699     1841
##    y1[3]    0.51   0.49 1.86 1.80  -2.54   3.48 1.00     1818     1954
##    y1[4]    4.98   5.07 1.84 1.82   1.92   7.90 1.00     1857     1986
##    y1[5]    1.56   1.54 1.88 1.86  -1.60   4.80 1.00     1878     1998
##    y1[6]    3.09   3.06 1.83 1.79   0.05   6.11 1.00     1892     1931
##    y1[7]    1.25   1.29 1.83 1.81  -1.79   4.23 1.00     1874     1797
##    y1[8]    1.24   1.26 1.88 1.88  -1.85   4.30 1.00     1978     2058
##    y1[9]    1.91   1.90 1.81 1.79  -0.96   4.92 1.00     1927     1701
##    y1[10]   4.64   4.68 1.92 1.91   1.52   7.73 1.00     1880     1808
##    y1[11]   1.04   1.05 1.80 1.76  -1.96   3.97 1.00     1866     1858
##    y1[12]   0.61   0.65 1.86 1.77  -2.54   3.65 1.00     2091     2030
##    y1[13]   2.70   2.71 1.82 1.79  -0.24   5.61 1.00     1562     1968
##    y1[14]   2.56   2.59 1.86 1.82  -0.44   5.60 1.00     1907     1873
##    y1[15]   1.51   1.50 1.82 1.74  -1.47   4.48 1.00     2001     1933
##    y1[16]   4.12   4.18 1.86 1.85   1.05   7.16 1.00     1924     1805
##    y1[17]  -1.23  -1.23 1.87 1.82  -4.33   1.80 1.00     1737     1768
##    y1[18]   0.49   0.51 1.80 1.73  -2.59   3.30 1.00     1906     1796
##    y1[19]   6.00   6.01 2.01 1.97   2.63   9.29 1.00     1919     1933
##    y1[20]   2.70   2.76 1.85 1.84  -0.45   5.70 1.00     1945     1776
##    y1[21]   1.80   1.85 1.85 1.80  -1.37   4.75 1.00     2026     1877
##    y1[22]   2.45   2.41 1.90 1.88  -0.60   5.55 1.00     1909     1846
##    y1[23]   2.17   2.15 1.87 1.89  -0.84   5.17 1.00     1909     1755
##    y1[24]  -0.45  -0.49 1.83 1.79  -3.43   2.54 1.00     1869     1972
##    y1[25]   0.59   0.59 1.90 1.89  -2.50   3.63 1.00     1852     1857
##    y1[26]   2.64   2.68 1.84 1.75  -0.56   5.64 1.00     1901     1846
##    y1[27]   3.83   3.84 1.90 1.89   0.73   6.92 1.00     1791     1900
##    y1[28]   2.79   2.84 1.86 1.90  -0.30   5.76 1.00     1758     1756
##    y1[29]  -0.21  -0.22 1.83 1.79  -3.27   2.84 1.00     1832     1973
##    y1[30]   0.65   0.65 1.85 1.80  -2.41   3.71 1.00     2010     2101
##    y1[31]   3.14   3.16 1.82 1.78   0.24   6.11 1.00     2091     1880
##    y1[32]   2.03   2.00 1.83 1.82  -1.04   4.99 1.00     1756     1956
##    y1[33]   1.49   1.56 1.81 1.72  -1.57   4.36 1.00     1650     1777
##    y1[34]   1.66   1.64 1.80 1.76  -1.34   4.60 1.00     1809     1969
##    y1[35]   3.54   3.56 1.91 1.88   0.40   6.70 1.00     1952     1890
##    y1[36]   2.95   2.97 1.90 1.88  -0.18   6.08 1.00     1915     2056
##    y1[37]   2.05   2.02 1.84 1.78  -1.05   5.01 1.00     2040     1900
##    y1[38]  -0.87  -0.88 1.80 1.75  -3.92   2.08 1.00     2019     1929
##    y1[39]   4.32   4.29 1.87 1.83   1.24   7.49 1.00     1982     2058
##    y1[40]   1.97   1.94 1.78 1.74  -0.89   4.89 1.00     1974     1956
##    y1[41]  -2.65  -2.62 1.93 1.91  -5.83   0.55 1.00     1719     2015
##    y1[42]   0.66   0.64 1.86 1.85  -2.39   3.75 1.00     1691     1596
##    y1[43]   2.90   2.94 1.84 1.83  -0.10   5.93 1.00     1921     1971
##    y1[44]  -1.89  -1.91 1.90 1.87  -4.96   1.25 1.00     2076     1917
##    y1[45]   2.74   2.72 1.83 1.72  -0.32   5.75 1.00     1738     1888
##    y1[46]   1.14   1.10 1.86 1.85  -1.80   4.24 1.00     1862     1928
##    y1[47]   1.41   1.43 1.84 1.79  -1.64   4.37 1.00     1885     1793
##    y1[48]  -0.48  -0.45 1.84 1.78  -3.47   2.44 1.00     2040     1732
##    y1[49]   3.69   3.60 1.97 1.97   0.54   6.94 1.00     1644     1800
##    y1[50]   1.60   1.55 1.82 1.84  -1.35   4.56 1.00     1862     1713
##    m1[1]    0.78   0.78 0.54 0.55  -0.10   1.64 1.00     1196     1199
##    m1[2]   -0.34  -0.33 0.45 0.44  -1.08   0.41 1.01      948      903
##    m1[3]    0.50   0.49 0.53 0.52  -0.38   1.38 1.00     1289     1187
##    m1[4]    5.02   5.02 0.63 0.62   4.00   6.06 1.00     1346     1086
##    m1[5]    1.56   1.56 0.53 0.51   0.71   2.42 1.00     1321     1287
##    m1[6]    3.05   3.05 0.51 0.50   2.23   3.87 1.00     1445     1356
##    m1[7]    1.25   1.25 0.57 0.57   0.33   2.18 1.00     1466     1404
##    m1[8]    1.31   1.31 0.65 0.63   0.18   2.36 1.00     1331     1375
##    m1[9]    1.90   1.90 0.53 0.53   1.03   2.76 1.00     1345     1453
##    m1[10]   4.63   4.65 0.64 0.63   3.60   5.67 1.00     1747     1476
##    m1[11]   1.04   1.02 0.49 0.49   0.24   1.84 1.00     1194     1317
##    m1[12]   0.56   0.56 0.58 0.58  -0.35   1.48 1.00     1689     1592
##    m1[13]   2.75   2.73 0.46 0.45   2.01   3.52 1.00     1184     1384
##    m1[14]   2.56   2.54 0.53 0.52   1.70   3.40 1.00     1303     1400
##    m1[15]   1.51   1.51 0.55 0.55   0.61   2.41 1.00     1272     1459
##    m1[16]   4.09   4.09 0.56 0.54   3.19   5.02 1.00     1338     1230
##    m1[17]  -1.21  -1.22 0.59 0.57  -2.21  -0.25 1.00     1228     1063
##    m1[18]   0.49   0.49 0.59 0.59  -0.44   1.43 1.00     1745     1568
##    m1[19]   5.98   5.96 0.80 0.79   4.68   7.28 1.00     1511     1284
##    m1[20]   2.68   2.69 0.59 0.58   1.73   3.60 1.00     1298     1052
##    m1[21]   1.73   1.73 0.61 0.62   0.71   2.71 1.00     1471     1451
##    m1[22]   2.43   2.42 0.61 0.63   1.45   3.42 1.01     1574     1450
##    m1[23]   2.25   2.26 0.51 0.50   1.43   3.07 1.00     1163     1131
##    m1[24]  -0.48  -0.49 0.49 0.48  -1.28   0.38 1.01     1121     1009
##    m1[25]   0.56   0.55 0.66 0.67  -0.49   1.64 1.00     1423     1441
##    m1[26]   2.62   2.61 0.50 0.51   1.81   3.45 1.00     1384     1485
##    m1[27]   3.84   3.83 0.55 0.54   2.94   4.75 1.00     1372     1299
##    m1[28]   2.87   2.86 0.46 0.45   2.13   3.62 1.00     1318     1454
##    m1[29]  -0.24  -0.25 0.50 0.49  -1.02   0.61 1.00     1348     1508
##    m1[30]   0.66   0.66 0.49 0.48  -0.14   1.47 1.00     1189     1224
##    m1[31]   3.14   3.13 0.55 0.56   2.25   4.02 1.01     2060     1424
##    m1[32]   1.99   2.00 0.51 0.49   1.15   2.80 1.00     1089     1193
##    m1[33]   1.53   1.53 0.43 0.43   0.82   2.23 1.01      966     1039
##    m1[34]   1.67   1.68 0.46 0.45   0.91   2.43 1.00     1013     1251
##    m1[35]   3.53   3.54 0.70 0.71   2.42   4.67 1.01     2423     1486
##    m1[36]   2.96   2.96 0.56 0.57   2.05   3.86 1.00     2024     1485
##    m1[37]   2.05   2.04 0.62 0.62   1.02   3.06 1.00     1884     1524
##    m1[38]  -0.85  -0.86 0.55 0.54  -1.74   0.09 1.00     1534     1485
##    m1[39]   4.43   4.45 0.52 0.52   3.57   5.27 1.00     1775     1560
##    m1[40]   1.95   1.96 0.51 0.49   1.10   2.78 1.00     1103     1194
##    m1[41]  -2.55  -2.57 0.78 0.77  -3.83  -1.27 1.00     1812     1554
##    m1[42]   0.68   0.68 0.47 0.48  -0.07   1.45 1.00     1143     1324
##    m1[43]   2.92   2.92 0.56 0.56   2.03   3.83 1.00     1513     1422
##    m1[44]  -1.88  -1.88 0.69 0.67  -3.03  -0.73 1.00     1696     1380
##    m1[45]   2.72   2.72 0.44 0.43   2.01   3.43 1.00     1427     1387
##    m1[46]   1.15   1.13 0.63 0.63   0.14   2.16 1.00     1500     1135
##    m1[47]   1.43   1.43 0.54 0.55   0.55   2.32 1.00     1419     1273
##    m1[48]  -0.43  -0.43 0.50 0.49  -1.28   0.36 1.00     1036      993
##    m1[49]   3.65   3.65 0.78 0.79   2.38   4.87 1.00     1253     1366
##    m1[50]   1.53   1.51 0.63 0.63   0.52   2.56 1.00     1542     1370

fn(f1)
## Initial log joint probability = -477.205 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       28      -46.7072   0.000155759    0.00087093           1           1       31    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -50.15 -49.80 2.10 1.98 -54.11 -47.46 1.00      688      964
##    b[1]     0.67   0.67 0.45 0.43  -0.06   1.40 1.00      714      810
##    b[2]     1.11   1.12 0.25 0.26   0.70   1.53 1.00     1891     1691
##    b[3]     1.08   1.08 0.25 0.23   0.66   1.48 1.00     1934     1280
##    b[4]     1.00   0.99 0.71 0.71  -0.19   2.17 1.01      505      652
##    b[5]     2.89   2.89 0.64 0.65   1.86   3.89 1.00      570      778
##    b[6]    -1.99  -1.94 1.03 1.03  -3.73  -0.38 1.00      411      527
##    s        1.70   1.68 0.19 0.18   1.42   2.04 1.00     1621     1556
##    y1[1]    1.33   1.31 1.80 1.82  -1.60   4.23 1.00     1716     1865
##    y1[2]   -0.84  -0.86 1.79 1.74  -3.75   2.24 1.00     1580     1930
##    y1[3]    0.09   0.08 1.82 1.87  -2.85   2.99 1.00     2004     1775
##    y1[4]    4.63   4.63 1.82 1.77   1.72   7.60 1.00     1900     1900
##    y1[5]    1.89   1.90 1.75 1.75  -0.98   4.82 1.00     1876     1786
##    y1[6]    3.49   3.51 1.75 1.73   0.63   6.27 1.00     1853     1969
##    y1[7]    0.55   0.57 1.81 1.85  -2.44   3.53 1.00     1831     1970
##    y1[8]    0.98   1.00 1.80 1.78  -1.96   3.94 1.00     1985     1857
##    y1[9]    2.28   2.26 1.81 1.75  -0.61   5.22 1.00     2018     1891
##    y1[10]   5.29   5.31 1.86 1.83   2.28   8.32 1.00     1805     1932
##    y1[11]   0.67   0.70 1.79 1.79  -2.31   3.58 1.00     1828     1685
##    y1[12]   1.09   1.11 1.79 1.80  -1.90   3.91 1.00     1729     1931
##    y1[13]   2.13   2.09 1.81 1.74  -0.87   5.13 1.00     2030     2015
##    y1[14]   1.93   1.96 1.80 1.72  -1.07   4.77 1.00     1896     1894
##    y1[15]   2.18   2.19 1.76 1.72  -0.65   5.05 1.00     1906     1604
##    y1[16]   3.54   3.52 1.79 1.68   0.59   6.54 1.00     2008     1785
##    y1[17]  -1.80  -1.76 1.78 1.81  -4.74   1.15 1.00     1772     2037
##    y1[18]   0.94   0.90 1.79 1.76  -2.00   3.93 1.00     1807     1921
##    y1[19]   5.65   5.68 1.96 1.91   2.28   8.82 1.00     1867     1957
##    y1[20]   2.45   2.43 1.78 1.79  -0.34   5.44 1.00     1819     1822
##    y1[21]   2.06   2.02 1.82 1.77  -0.87   5.04 1.00     1807     1973
##    y1[22]   2.79   2.86 1.80 1.81  -0.22   5.67 1.00     2054     1697
##    y1[23]   2.01   1.99 1.76 1.77  -0.86   4.88 1.00     1970     1762
##    y1[24]  -0.91  -0.88 1.77 1.75  -3.77   1.95 1.00     1625     1855
##    y1[25]   1.11   1.14 1.85 1.81  -2.05   4.13 1.00     1728     1917
##    y1[26]   3.06   3.06 1.76 1.79   0.14   5.87 1.00     1841     2036
##    y1[27]   3.26   3.30 1.83 1.83   0.30   6.20 1.00     1979     1856
##    y1[28]   3.37   3.34 1.78 1.65   0.45   6.31 1.00     1809     1856
##    y1[29]   0.14   0.16 1.79 1.76  -2.89   3.09 1.00     1805     1889
##    y1[30]   0.29   0.29 1.74 1.69  -2.58   3.12 1.00     1962     1891
##    y1[31]   3.52   3.54 1.78 1.70   0.52   6.43 1.00     1940     1930
##    y1[32]   1.71   1.67 1.76 1.75  -1.25   4.58 1.00     1909     1973
##    y1[33]   1.19   1.22 1.75 1.74  -1.70   3.99 1.00     1985     1981
##    y1[34]   1.31   1.30 1.71 1.63  -1.45   4.14 1.00     1720     2105
##    y1[35]   4.07   4.08 1.85 1.87   1.08   7.06 1.00     2037     1973
##    y1[36]   3.46   3.42 1.86 1.86   0.43   6.55 1.00     1502     1420
##    y1[37]   2.39   2.37 1.87 1.80  -0.67   5.47 1.00     1802     1930
##    y1[38]  -0.39  -0.39 1.80 1.80  -3.32   2.46 1.00     1706     1706
##    y1[39]   4.98   4.97 1.78 1.69   2.06   7.98 1.00     1908     1883
##    y1[40]   1.64   1.60 1.78 1.74  -1.30   4.56 1.00     1794     1924
##    y1[41]  -2.12  -2.08 1.91 1.93  -5.25   1.10 1.00     1813     1823
##    y1[42]   1.24   1.23 1.80 1.70  -1.78   4.25 1.00     1878     1721
##    y1[43]   2.27   2.32 1.79 1.78  -0.65   5.27 1.00     1850     1839
##    y1[44]  -1.55  -1.55 1.76 1.74  -4.41   1.40 1.00     2098     1843
##    y1[45]   3.13   3.14 1.79 1.77   0.14   6.10 1.00     1972     1893
##    y1[46]   0.38   0.42 1.84 1.90  -2.69   3.39 1.00     1728     1840
##    y1[47]   0.79   0.78 1.82 1.83  -2.20   3.88 1.00     2096     1904
##    y1[48]  -0.92  -0.97 1.80 1.78  -3.80   2.06 1.00     1711     1827
##    y1[49]   4.41   4.46 1.94 1.91   1.29   7.68 1.00     1624     1882
##    y1[50]   1.21   1.22 1.80 1.80  -1.72   4.12 1.00     1924     1815
##    m1[1]    1.36   1.35 0.60 0.59   0.35   2.37 1.00      898     1187
##    m1[2]   -0.81  -0.82 0.51 0.49  -1.64   0.05 1.00      661      947
##    m1[3]    0.10   0.09 0.56 0.56  -0.80   1.01 1.00      934     1314
##    m1[4]    4.59   4.60 0.66 0.64   3.50   5.65 1.00     1925     1486
##    m1[5]    1.91   1.91 0.54 0.52   1.02   2.76 1.00     1668     1592
##    m1[6]    3.50   3.50 0.52 0.51   2.63   4.33 1.00     1065     1337
##    m1[7]    0.53   0.54 0.66 0.64  -0.55   1.62 1.00     1080     1138
##    m1[8]    0.96   0.94 0.66 0.66  -0.10   2.03 1.00     1200     1230
##    m1[9]    2.26   2.27 0.54 0.54   1.38   3.14 1.00     1294     1543
##    m1[10]   5.20   5.21 0.67 0.66   4.11   6.32 1.00     1192     1463
##    m1[11]   0.68   0.68 0.52 0.51  -0.16   1.53 1.00      904     1191
##    m1[12]   1.11   1.10 0.60 0.59   0.12   2.12 1.00     1192     1445
##    m1[13]   2.15   2.15 0.55 0.53   1.25   3.05 1.01     1131     1439
##    m1[14]   1.94   1.93 0.60 0.59   0.95   2.91 1.00     1196     1444
##    m1[15]   2.13   2.14 0.61 0.61   1.12   3.17 1.00      931     1177
##    m1[16]   3.59   3.59 0.60 0.58   2.58   4.54 1.00     1536     1354
##    m1[17]  -1.75  -1.76 0.65 0.61  -2.83  -0.66 1.00      771     1133
##    m1[18]   1.04   1.03 0.61 0.61   0.03   2.07 1.00     1242     1462
##    m1[19]   5.62   5.63 0.80 0.77   4.32   6.93 1.01     2349     1624
##    m1[20]   2.45   2.45 0.58 0.57   1.49   3.45 1.00     1263     1218
##    m1[21]   2.07   2.07 0.60 0.62   1.08   3.05 1.00     1430     1573
##    m1[22]   2.83   2.83 0.61 0.60   1.81   3.81 1.00     1335     1673
##    m1[23]   1.98   1.98 0.52 0.50   1.14   2.86 1.00     1085     1137
##    m1[24]  -0.95  -0.96 0.54 0.53  -1.82  -0.02 1.00      730     1097
##    m1[25]   1.12   1.11 0.70 0.69  -0.03   2.29 1.00     1099     1380
##    m1[26]   3.04   3.04 0.52 0.51   2.19   3.87 1.00     1104     1388
##    m1[27]   3.31   3.33 0.60 0.58   2.30   4.28 1.00     1485     1222
##    m1[28]   3.30   3.30 0.48 0.46   2.50   4.06 1.00     1002     1182
##    m1[29]   0.25   0.25 0.53 0.51  -0.63   1.11 1.00     1139     1287
##    m1[30]   0.28   0.28 0.52 0.52  -0.57   1.14 1.00      863     1228
##    m1[31]   3.60   3.58 0.58 0.59   2.67   4.57 1.00     1472     1313
##    m1[32]   1.69   1.69 0.52 0.52   0.87   2.58 1.01     1046     1298
##    m1[33]   1.20   1.20 0.46 0.46   0.45   1.95 1.00      821     1036
##    m1[34]   1.36   1.35 0.49 0.49   0.59   2.17 1.01      901     1085
##    m1[35]   4.03   4.02 0.72 0.71   2.86   5.22 1.00     1679     1512
##    m1[36]   3.41   3.39 0.58 0.59   2.47   4.38 1.00     1536     1352
##    m1[37]   2.43   2.42 0.63 0.64   1.38   3.46 1.00     1816     1562
##    m1[38]  -0.40  -0.41 0.56 0.55  -1.33   0.50 1.00     1398     1323
##    m1[39]   4.99   4.98 0.57 0.57   4.04   5.93 1.00     1036     1222
##    m1[40]   1.66   1.66 0.53 0.52   0.82   2.56 1.01     1051     1335
##    m1[41]  -2.24  -2.25 0.75 0.76  -3.48  -1.03 1.00     2032     1261
##    m1[42]   1.24   1.23 0.53 0.51   0.39   2.12 1.00      865     1047
##    m1[43]   2.32   2.34 0.62 0.60   1.26   3.31 1.00     1318     1133
##    m1[44]  -1.51  -1.51 0.69 0.69  -2.66  -0.40 1.00     1711     1303
##    m1[45]   3.15   3.14 0.47 0.47   2.38   3.94 1.00     1130     1402
##    m1[46]   0.42   0.43 0.71 0.70  -0.71   1.58 1.00     1185     1434
##    m1[47]   0.72   0.73 0.64 0.62  -0.32   1.79 1.00     1076     1130
##    m1[48]  -0.91  -0.91 0.55 0.53  -1.82   0.01 1.00      715     1128
##    m1[49]   4.44   4.44 0.86 0.87   3.04   5.86 1.00      902     1073
##    m1[50]   1.21   1.21 0.63 0.63   0.18   2.24 1.00     1209     1555

fn(f2)
## Initial log joint probability = -123.462 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       35       -22.216   6.34509e-05    0.00153826      0.2246           1       42    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -26.78 -26.42 2.24 2.12 -30.98 -23.80 1.00      773     1160
##    b[1]     0.23   0.23 0.28 0.27  -0.23   0.70 1.00      980     1464
##    b[2]     0.96   0.95 0.15 0.15   0.71   1.21 1.00     2509     1773
##    b[3]     0.84   0.84 0.15 0.14   0.59   1.09 1.00     2380     1292
##    b[4]     1.46   1.45 0.44 0.44   0.75   2.18 1.00      756      953
##    b[5]     2.20   2.20 0.41 0.41   1.55   2.87 1.00      618      914
##    b[6]    -1.15  -1.15 0.14 0.14  -1.38  -0.91 1.00     2011     1770
##    b[7]    -1.13  -1.14 0.65 0.68  -2.23  -0.04 1.00      483      612
##    s        1.05   1.04 0.12 0.12   0.87   1.27 1.00     2022     1496
##    y1[1]    1.90   1.91 1.14 1.12   0.03   3.77 1.00     1962     2014
##    y1[2]   -1.49  -1.46 1.12 1.08  -3.35   0.31 1.00     1924     1931
##    y1[3]    1.15   1.14 1.09 1.10  -0.64   2.97 1.00     1902     2012
##    y1[4]    3.45   3.46 1.16 1.18   1.52   5.33 1.00     1816     1894
##    y1[5]    0.65   0.66 1.15 1.16  -1.27   2.50 1.00     2010     2058
##    y1[6]    3.35   3.36 1.11 1.11   1.53   5.16 1.00     1967     1997
##    y1[7]    0.25   0.24 1.12 1.10  -1.57   2.04 1.00     1628     1792
##    y1[8]    3.54   3.56 1.19 1.20   1.63   5.45 1.00     1974     1846
##    y1[9]    1.73   1.71 1.11 1.12  -0.15   3.53 1.00     1964     1881
##    y1[10]   4.69   4.67 1.14 1.09   2.82   6.59 1.00     1889     1911
##    y1[11]   1.12   1.10 1.10 1.11  -0.76   2.97 1.00     1563     1845
##    y1[12]   3.06   3.07 1.16 1.09   1.07   5.06 1.00     2131     2012
##    y1[13]   2.39   2.40 1.09 1.08   0.64   4.17 1.00     1927     1974
##    y1[14]   2.60   2.56 1.16 1.10   0.73   4.52 1.00     1974     1971
##    y1[15]   2.95   2.94 1.14 1.09   1.07   4.78 1.00     1569     1495
##    y1[16]   4.04   4.03 1.12 1.05   2.15   5.91 1.00     1810     1800
##    y1[17]  -2.54  -2.55 1.13 1.10  -4.39  -0.69 1.00     1735     1894
##    y1[18]   3.14   3.15 1.18 1.16   1.19   5.06 1.00     2034     1885
##    y1[19]   3.37   3.34 1.18 1.16   1.43   5.38 1.00     2286     2056
##    y1[20]   1.27   1.27 1.11 1.07  -0.52   3.06 1.00     2125     1872
##    y1[21]   2.32   2.31 1.17 1.15   0.38   4.21 1.00     1827     1871
##    y1[22]   3.95   3.98 1.15 1.14   2.01   5.79 1.00     2097     1660
##    y1[23]   0.90   0.88 1.08 1.06  -0.88   2.66 1.00     1920     2046
##    y1[24]  -1.41  -1.42 1.13 1.15  -3.25   0.44 1.00     1778     1955
##    y1[25]   3.17   3.20 1.15 1.13   1.26   5.07 1.00     1733     1626
##    y1[26]   2.80   2.77 1.10 1.06   1.07   4.62 1.00     1885     1700
##    y1[27]   4.18   4.16 1.12 1.11   2.36   6.04 1.00     1979     2030
##    y1[28]   2.54   2.56 1.09 1.06   0.73   4.31 1.00     1868     1477
##    y1[29]   0.09   0.11 1.07 1.01  -1.69   1.87 1.00     1978     2054
##    y1[30]   0.76   0.74 1.10 1.04  -1.03   2.61 1.00     1971     1937
##    y1[31]   4.05   4.01 1.13 1.13   2.16   5.94 1.00     1967     1973
##    y1[32]   1.23   1.22 1.12 1.11  -0.56   3.07 1.00     1832     1757
##    y1[33]   0.63   0.61 1.11 1.09  -1.21   2.51 1.00     1978     1971
##    y1[34]   0.89   0.89 1.10 1.06  -0.91   2.71 1.00     1966     1876
##    y1[35]   6.46   6.44 1.18 1.13   4.51   8.36 1.00     2022     1908
##    y1[36]   3.98   3.97 1.09 1.04   2.19   5.74 1.00     2018     1970
##    y1[37]   3.36   3.36 1.12 1.05   1.44   5.20 1.00     1953     2011
##    y1[38]  -0.91  -0.93 1.11 1.09  -2.72   0.89 1.00     2048     2035
##    y1[39]   3.35   3.37 1.15 1.10   1.40   5.22 1.00     1715     1911
##    y1[40]   1.30   1.29 1.11 1.10  -0.46   3.09 1.00     1720     1910
##    y1[41]  -5.17  -5.19 1.22 1.18  -7.12  -3.16 1.00     1987     1967
##    y1[42]   1.36   1.37 1.11 1.11  -0.43   3.14 1.00     1824     2014
##    y1[43]   4.13   4.11 1.17 1.20   2.22   6.08 1.00     1784     1835
##    y1[44]  -3.24  -3.24 1.16 1.17  -5.13  -1.34 1.00     2029     1981
##    y1[45]   2.25   2.28 1.12 1.07   0.37   4.06 1.00     1755     1704
##    y1[46]   0.40   0.41 1.16 1.11  -1.48   2.30 1.00     2091     1849
##    y1[47]   0.53   0.51 1.15 1.14  -1.35   2.37 1.00     1679     1787
##    y1[48]  -1.15  -1.13 1.09 1.09  -2.97   0.57 1.00     1862     2051
##    y1[49]   2.11   2.12 1.18 1.18   0.22   4.02 1.00     1996     1827
##    y1[50]   3.12   3.11 1.15 1.12   1.19   5.00 1.00     1659     1887
##    m1[1]    1.89   1.88 0.36 0.38   1.30   2.48 1.00     1334     1392
##    m1[2]   -1.51  -1.50 0.33 0.32  -2.05  -0.95 1.00     1048     1497
##    m1[3]    1.12   1.11 0.37 0.36   0.53   1.74 1.00     1073     1495
##    m1[4]    3.44   3.44 0.43 0.43   2.75   4.15 1.00     2613     1607
##    m1[5]    0.63   0.64 0.36 0.36   0.05   1.22 1.00     1272     1439
##    m1[6]    3.34   3.34 0.34 0.33   2.78   3.90 1.00     1366     1567
##    m1[7]    0.24   0.24 0.41 0.41  -0.42   0.93 1.00     1337     1595
##    m1[8]    3.53   3.52 0.53 0.50   2.66   4.40 1.00     1435     1341
##    m1[9]    1.69   1.69 0.35 0.34   1.11   2.26 1.00     1392     1460
##    m1[10]   4.72   4.73 0.43 0.42   4.02   5.41 1.00     1397     1531
##    m1[11]   1.13   1.13 0.32 0.31   0.62   1.68 1.00     1035     1411
##    m1[12]   3.05   3.05 0.45 0.42   2.32   3.79 1.00     1800     1579
##    m1[13]   2.37   2.37 0.35 0.35   1.79   2.93 1.00     1355     1407
##    m1[14]   2.63   2.62 0.39 0.40   2.01   3.25 1.00     1380     1422
##    m1[15]   2.97   2.96 0.39 0.39   2.34   3.60 1.00     1376     1394
##    m1[16]   4.03   4.03 0.39 0.38   3.39   4.67 1.00     1916     1528
##    m1[17]  -2.56  -2.56 0.42 0.42  -3.27  -1.88 1.00     1240     1350
##    m1[18]   3.18   3.17 0.47 0.44   2.42   3.95 1.00     1861     1603
##    m1[19]   3.33   3.33 0.57 0.57   2.42   4.26 1.00     2865     1670
##    m1[20]   1.26   1.26 0.37 0.36   0.65   1.87 1.00     1808     1207
##    m1[21]   2.31   2.31 0.38 0.38   1.68   2.93 1.00     1936     1627
##    m1[22]   3.96   3.97 0.41 0.41   3.27   4.61 1.00     2410     1676
##    m1[23]   0.91   0.91 0.33 0.33   0.38   1.45 1.00     1458     1238
##    m1[24]  -1.44  -1.44 0.34 0.34  -2.01  -0.87 1.00     1107     1584
##    m1[25]   3.14   3.13 0.49 0.49   2.36   3.96 1.00     1693     1617
##    m1[26]   2.80   2.81 0.33 0.33   2.26   3.34 1.00     1370     1522
##    m1[27]   4.20   4.20 0.40 0.40   3.53   4.86 1.00     1727     1619
##    m1[28]   2.58   2.57 0.32 0.31   2.04   3.12 1.00     1125     1349
##    m1[29]   0.11   0.11 0.33 0.32  -0.43   0.66 1.00     1431     1488
##    m1[30]   0.81   0.80 0.33 0.32   0.29   1.36 1.00     1026     1370
##    m1[31]   4.07   4.06 0.36 0.35   3.48   4.67 1.00     1688     1585
##    m1[32]   1.23   1.23 0.32 0.33   0.71   1.75 1.00     1239     1375
##    m1[33]   0.61   0.61 0.29 0.28   0.14   1.08 1.00     1095     1297
##    m1[34]   0.90   0.90 0.30 0.30   0.41   1.39 1.00     1122     1354
##    m1[35]   6.47   6.48 0.54 0.53   5.60   7.37 1.00     2443     1579
##    m1[36]   3.99   3.98 0.36 0.36   3.39   4.60 1.00     1762     1585
##    m1[37]   3.35   3.36 0.40 0.39   2.69   3.99 1.00     2167     1536
##    m1[38]  -0.90  -0.90 0.35 0.35  -1.48  -0.32 1.00     1678     1623
##    m1[39]   3.38   3.38 0.41 0.41   2.70   4.03 1.00     1155     1422
##    m1[40]   1.31   1.30 0.32 0.32   0.78   1.84 1.00     1220     1388
##    m1[41]  -5.15  -5.13 0.59 0.59  -6.13  -4.20 1.00     2175     1356
##    m1[42]   1.40   1.39 0.33 0.33   0.88   1.95 1.00     1236     1241
##    m1[43]   4.13   4.13 0.46 0.46   3.39   4.89 1.00     1469     1255
##    m1[44]  -3.21  -3.21 0.47 0.48  -4.00  -2.47 1.00     2061     1389
##    m1[45]   2.20   2.20 0.31 0.31   1.70   2.71 1.00     1082     1375
##    m1[46]   0.37   0.36 0.43 0.43  -0.32   1.10 1.00     1391     1530
##    m1[47]   0.54   0.53 0.39 0.39  -0.10   1.20 1.00     1314     1518
##    m1[48]  -1.14  -1.14 0.35 0.34  -1.71  -0.56 1.00     1052     1222
##    m1[49]   2.15   2.14 0.59 0.58   1.17   3.13 1.00     1364     1445
##    m1[50]   3.12   3.12 0.46 0.45   2.37   3.90 1.00     1399     1443